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A holding  time  model  of  serial  production  lines  with  an 
arbitrary  number  of  stations  and  an  arbitrary  distribution 
of  service  times  is  presented.  A system  of  integral  equa- 
tions for  the  distribution  of  interdeparture  time  is 
derived.  For  a special  class  of  distributions  of  service 
times,  two  numerical  procedures,  one  exact  and  the  other 
approximate,  for  calculating  the  line  throughput  rate  are 
given.  For  the  three-station  case  the  exact  procedure  is 
specialized  into  a formula  for  the  throughput  rate  of  a 
three-station  production  line.  This  formula  is  used  to 
investigate  the  effect  of  line  unbalance  on  the  throughput 
rate.  A model  of  production  lines  with  station  breakdown 
and  repair  is  also  considered  and  a method  to  obtain  an 
equivalent  production  line  without  station  breakdown  is 
derived . 
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CHAPTER  1 


INTRODUCTION 


1.1.  General  Description  of  a Production  Line 

The  production  line  is  a common  material  handling  and 
processing  device  of  modern  manufacturing  systems,  and  the 
design  of  efficient  production  lines  is  one  of  the  most 
important  jobs  of  industrial  engineers. 

The  production  line  model  to  be  considered  here  is  a 
series  arrangement  of  work  stations.  At  each  station  a 
specific  operation  is  performed  on  the  workpiece  passing 
through  the  station.  Each  workpiece  passes  through  all  of 
the  stations  in  the  same  sequence.  An  important  measure  of 
the  line's  efficiency  is  its  throughput  rate  or  mean  pro- 
duction rate  which  is  defined  to  be  the  number  of  workpieces 
released  from  the  line  per  unit  time  in  the  long  run. 

An  ideal  production  line  is  one  in  which 

- the  station  service  time  is  deterministic 

- the  line  is  balanced,  which  means  that  every  station 
processes  the  same  number  of  items  per  unit  time 

- there  is  no  station  breakdown. 

Actual  production  lines  suffer  losses  in  efficiency 
owing  to  three  distinct  causes.  These  causes  are 
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- randomness  of  the  station  service  times 

- unbalance  of  the  line 

- station  breakdown. 

Each  of  the  causes  results  in  a reduction  in  the  line 
efficiency.  Losses  in  efficiency,  as  observed  on  a par- 
ticular station,  are  evidenced  by  periods  where  that  station 
is  either  broken  down,  blocked  or  starved.  A station  may 
break  down  and  cease  to  operate  because  of  a direct  station 
failure,  or  it  may  have  to  be  shut  down  because  the  perfor- 
mance of  the  station  is  not  within  the  bounds  of  prescribed 
tolerances.  A station  is  blocked  if  it  has  completed 
service  but  cannot  pass  on  its  workpiece  to  the  next  station 
because  that  station  is  busy,  blocked  or  broken  down.  A 
station  is  starved  if  it  has  passed  on  its  workpiece  to  the 
next  station  but  has  not  yet  received  a new  workpiece  from 
the  preceding  station.  Blocking  and  starving  are  caused  by 
unbalance,  randomness  of  station  service  times  and  station 
breakdown.  The  effect  of  blocking  propagates  backward 
through  the  line,  while  the  effect  of  starving  propagates 
forward.  A station  breakdown  will  cause  preceding  stations 
to  be  blocked  and  succeeding  stations  to  be  starved.  Thus 
there  is  a strong  coupling  effect  between  stations. 

The  coupling  effect  can  be  reduced  by  placing  buffers, 
in-process  storage,  between  stations.  Buffers  have  a 
decoupling  effect  since  an  otherwise  starved  station  can 
still  be  fed  from  the  preceding  buffer,  as  long  as  this 
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buffer  is  not  empty.  Similarly,  a blocked  station  can 
release  its  finished  workpiece  into  the  subsequent  buffer, 
as  long  as  this  buffer  is  not  full.  A finite  number  of  buf- 
fer does  not  completely  decouple  adjacent  stations.  Station 
interference  still  exists  when  the  buffer  is  either  empty  or 
full.  In  practice  buffers  are  of  small  size.  A buffer  with 
a capacity  for  L items  is  equivalent  to  L stations  with  zero 
service  times  in  series.  Therefore,  in  principle,  the  case 
of  a line  with  buffers  does  not  need  to  be  treated  differ- 
ently from  that  of  a line  without  buffers. 

The  throughput  rate  is  a function  of 

- the  number  of  stations 

- the  distribution  of  station  service  times 

- the  distribution  of  station  repair  times 

- the  distribution  of  station  uptime  between  breakdowns 
In  this  dissertation  we  investigate  models  that  incorporate 
these  functional  relationships  and  we  develop  methods  for 
computing  the  throughput  rate. 

1.2.  Literature  Review 

Research  on  production  lines  appeared  in  the  literature 
in  the  early  1950s.  According  to  Buzacott  [1972]  , produc- 
tion lines  were  first  studied  analytically  via  a probabilis- 
tic approach  by  Vladzievskii  [1952] . 

One  common  probabilistic  model  for  studying  production 
lines  is  the  system  state  model.  The  various  states  that 
the  system  may  assume  and  the  transition  characteristics 
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among  them  are  identified.  Determination  of  throughput  rate 
is  based  on  finding  the  steady  state  probabilities  for  those 
states  in  which  the  last  station  is  actively  processing 
workpieces . 

Simulation  models  have  been  applied  with  the  expecta- 
tion of  obtaining  some  of  the  properties  of  certain  produc- 
tion lines.  It  seems  unlikely  though  that  simulation  will 
uncover  general  relationships.  Nevertheless,  simulation  can 
be  considered  as  a useful  tool  for  verifying  and  comparing 
different  models. 

Continuous  models  have  also  been  used  to  represent 
production  lines.  In  these  models  workpieces  are  treated  as 
a continuous  fluid  rather  than  as  discrete  items. 

Individual  studies  differ  by  the  model  used  to 
represent  a production  line  and  by  the  method  of  analysis. 
One  major  difference  is  whether  or  not  station  breakdown  is 
considered.  The  case  of  station  breakdown  is  further 
differentiated  by  considering  the  station  service  time  to  be 
either  random  or  fixed.  We  divide  the  literature  into  three 
major  classes,  according  to  the  model  used.  The  first  class 
consists  of  papers  that  treat  models  of  production  lines  in 
which  station  breakdown  is  not  considered.  The  second  class 
consists  of  papers  that  treat  models  of  production  lines  in 
which  stations  break  down  at  random  and  station  service 
times  are  fixed  and  identical.  The  third  class  consists  of 
papers  that  treat  models  of  production  lines  in  which 
stations  break  down  at  random  and  station  service  times  are 
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random  variables.  The  literature  review  that  follows  is 
divided  into  these  three  classes.  We  only  discuss  papers 
that  use  analytical  methods  that  are  based  on  discrete 
models  of  production  lines,  since  these  are  related  to  the 
model  and  methods  of  analysis  of  this  dissertation.  Simu- 
lation models  can  be  found  in  Freeman  [1964]  , Anderson  and 
Moodie  [1969],  Masso  and  Smith  [1974]  and  Ho  et  al.  [1979]. 
Continuous  models  can  be  found  in  Sevast'yanov  [1962]  , 

Murphy  [1978]  , Wijngaard  [1979]  and  Gershwin  and  Schick 
[1980]  . 

Unless  otherwise  stated,  all  papers  reviewed  here 
assume  that  there  is  an  infinite  supply  of  workpieces  to  the 
first  station  and  that  there  is  always  space  available  into 
which  the  last  station  can  discharge  a finished  item. 

1.2.1.  Models  of  Production  Lines  without  Station  Breakdown 

Hunt  [1956]  considers  production  lines  with  exponential 
service  times.  Using  a Markov  model  he  derives  a simple 
closed  form  expression  for  the  throughput  rate  of  a two- 
station  balanced  production  line  with  an  arbitrary  number  of 
buffers  between  stations.  He  also  derives  a closed  form 
expression  for  the  throughput  rate  of  a three-station 
production  line  with  no  buffers.  Even  though  the  throughput 
^3te  is  a function  of  only  three  variables,  the  expression 
turns  out  to  be  very  complicated.  It  has  a numerator 
polynomial  of  degree  8 with  22  terms  and  a denominator 
polynomial  of  degree  7 with  24  terms.  This  complicated 
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formula  is  an  indication  that  it  is  practically  impossible 
to  obtain  closed  form  expressions  for  the  throughput  rate  of 
a production  line  with  more  than  three  stations.  Hunt  also 
gives  a numerical  solution  for  a four-station  production 
line  with  no  buffers  and  a three-station  balanced  production 
line  with  buffers  of  capacity  one. 

Avi-Itzak  and  Yadin  [1965]  study  a two-station  produc- 
tion line  with  general  service  times  and  with  no  buffers. 
Items  arrive  at  this  production  line  according  to  a Poisson 
process.  The  authors  obtain  the  moment  generating  function 
for  the  random  variable  time  spent  in  the  line  by  an  item. 

Hillier  and  Boling  [1966] , using  the  result  of  Hunt 
[1956] , show  numerically  that  the  throughput  rate  of  a 
three-station  production  line  with  exponential  service  time 
can  be  increased  by  assigning  a lower  mean  service  times  to 
the  middle  station  than  to  the  two  outside  stations.  They 
refer  to  this  allocation  of  mean  service  times  as  "bowl 
phenomenon . " 

Hillier  and  Boling  [1967],  using  a Markov  model, 
develop  a numerical  procedure  for  calculating  the  line 
throughput  rate.  They  treat  lines  in  which  all  stations  are 
identical  and  have  the  same  number  of  buffers,  for  the  case 
of  exponential  and  Erlang  service  times.  Their  solution 
method  exploits  the  sparsity  of  the  associated  transition 
rate  matrix.  Their  numerical  results  are  limited  to  at  most 
six  stations  and  at  most  four  buffers  per  station.  They 
also  give  an  approximation  method  for  calculating  the 


throughput  rate  of  a production  line  with  exponential 
service  times.  This  method  is  based  on  Burke's  theorem. 
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Muth  [1973]  considers  a production  line  with  an  arbi- 
trary distribution  of  service  times  and  an  arbitrary  number 
of  stations.  He  gives  methods  for  calculating  an  upper  and 
a lower  bounds  on  the  throughput  rate.  He  also  derives  an 
integral  equation  for  the  distribution  of  the  interdemand 
times  of  a three-station  production  line. 

Rao  [1975a]  uses  the  model  introduced  by  Muth  [1973]  to 
calculate  the  throughput  rate  of  two-station  production 
lines  with  Erlang,  uniform  and  normal  service  times.  In  a 
second  paper,  Rao  [1975b]  uses  his  previous  result  to 
develop  a method  for  calculating  throughput  rates  of  two- 
station  production  lines  with  an  arbitrary  number  of  buf- 
fers, where  the  first  station  has  exponential  service  times 
and  the  second  station  has  Erlang,  uniform  or  normal  service 
times . 

Using  the  integral  formulation  for  a three-station 
production  line  given  by  Muth  [1973]  , Rao  [1976]  derives 
closed  form  expressions  for  throughput  rates  of  several 
models  of  production  lines.  In  the  first  model  two  stations 
have  exponential  service  times  and  one  has  a deterministic 
service  time,  in  the  second  model  two  stations  have  deter- 
ministic service  times  and  one  has  an  exponential  service 
time,  and  in  the  third  model  the  two  outside  stations  have 
exponential  service  times  and  the  middle  station  has  a 
uniform  service  time.  Using  these  closed  form  expressions. 
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he  shows  numerically  that  the  bowl  phenomenon  does  not 
always  occur  in  a production  line  where  stations  have  widely 
different  coefficients  of  variation  of  service  times.  He 
finds  that  the  throughput  rate  can  be  increased  by  assigning 
a lower  mean  service  times  to  stations  with  higher  varia- 
bilities. He  calls  this  allocation  "variability  imbalance." 
He  concludes  that  for  some  cases  the  effect  of  variability 
imbalance  is  opposite  to  the  effect  of  mean  value  unbalance 
(bowl  phenomenon  of  Hillier  and  Boling  [1966])  and  that  the 
two  effects  may  cancel  one  another. 

Muth  [1977],  as  in  his  earlier  paper,  considers  a 
production  line  with  an  arbitrary  distribution  of  service 
times  and  an  arbitrary  number  of  stations.  He  uses  a 
passage  time  model  to  derive  relations  among  the  random 
variables  involved.  For  the  three-station  case  he  derives  a 
system  of  integral  equations  for  the  distribution  of  the 
interdemand  times.  He  also  gives  two  numerical  procedures 
for  solving  this  system  of  integral  equations  in  the  case  of 
identical  uniform  service  times  and  identical  Erlang  service 
times . 

Hillier  and  Boling  [1979],  using  the  numerical  pro- 
cedure they  developed  in  their  earlier  paper  [1967] , 
describe  the  bowl  phenomenon  for  symmetrical  production 
lines  with  up  to  six  stations  and  coefficient  of  variation 
of  service  times  equal  to  l//in  where  m=l,  2,  . . .,  7. 

They  assume  that  the  service  times  of  all  stations  are 
Erlang  with  the  same  coefficient  of  variation.  They  show 
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that  in  all  cases  considered,  the  line  throughput  rate  can 
be  increased  by  assigning  mean  service  times  that  decrease 
toward  the  middle  of  the  line. 

Altiok  [1982]  proposes  an  approximation  method  for 
calculating  the  throughput  rate  of  a production  line  where 
all  stations  have  exponential  service  times.  His  method  is 
based  on  the  no-memory  property  of  exponential  random 
variables  and  the  assumption  that  two  consecutive  stations 
are  never  blocked  at  the  same  time.  His  method,  as  with 
Hillier  and  Boling's  approximation  method,  gives  a good 
result  only  when  there  are  buffers  of  large  size  placed 
between  stations. 

Muth  [1984]  uses  a passage  time  model  to  give  activity 
network  representations  of  a production  line  with  an  arbi- 
trary distribution  of  service  times  and  an  arbitrary  number 
of  stations.  These  networks  are  a useful  tool  for  deriving 
relationships  among  random  variable  involved.  He  gives  a 
system  of  integral  equations  for  the  distribution  of  the 
interdemand  times  of  a three-station  production  line.  In 
the  case  all  stations  have  exponential  service  times,  these 
equations  are  reduced  into  a simple  closed  form  expression 
compared  to  Hunt's.  He  also  derives  a numerical  procedure 
for  the  case  of  two-stage  general  Erlang  service  times. 

1.2.2.  Models  of  Production  Lines  with  Station  Breakdown 
and  Fixed  Service  Times  ~~  ~ 

Buzacott  [1967]  derives  an  exact  formula  for  the 
throughput  rate  of  a two-station  production  line  with  finite 
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buffers.  For  a three-station  line  with  identical  repair 
time,  an  iterative  solution  procedure  for  finding  an  approx- 
imate throughput  rate  is  given.  The  extension  of  this 
iterative  procedure  to  lines  with  more  than  three  stations 
is  discussed.  The  iterative  approximation  is  based  on  a 
process  of  dividing  the  production  line  into  two  stages  and 
applying  the  formula  derived  for  the  two-station  line.  He 
also  gives  some  insight  on  how  and  where  to  divide  the  line 
into  stages.  In  a second  paper,  Buzacott  [1968]  considers 
several  line  configurations  where  stations  are  arranged  in 
any  combination  of  parallel-series  and  series-parallel  forms 
and  where  there  are  either  no  buffers  or  there  is  an  infi- 
nite number  of  buffers. 

Sheskin  [1976]  assumes  geometrically  distributed  repair 
times  and  time  between  breakdowns.  These  assumptions  give 
rise  to  a discrete  parameter  Markov  process.  He  also 
assumes  that  starved  or  blocked  stations  can  break  down  and 
= 1/  where  is  the  conditional  probability  that 
station  i breaks  down  at  cycle  n + 1 given  that  it  was 
operational  during  cycle  n,  and  r^  is  the  conditional 
probability  that  station  i is  repaired  during  cycle  n + 1 
given  that  it  was  under  repair  during  cycle  n.  Exact 
numerical  solutions  for  the  steady  state  probabilities  of 
the  resulting  Markov  chain  are  obtained  by  a method  he  calls 
"compression"  of  the  transition  probability  matrix. 
Approximate  numerical  solutions  are  obtained  by  a method  he 
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calls  "decomposition."  Using  these  methods,  he  is  able  to 
calculate  the  throughput  rate  of  a four-station  production 
line  containing  up  to  twelve  buffers. 

Okamura  and  Yamashima  [1977]  assume  geometrically  dis- 
tributed repair  times  and  time  between  breakdowns.  They 
compute  the  throughput  rate  of  a two-station  production  line 
by  solving  a system  of  (4M  + 6)  simultaneous  equations, 
where  M is  the  number  of  buffers.  Their  numerical  pro- 
cedure, unlike  the  approach  of  Hillier  and  Boling  [1967], 
does  not  exploit  the  structure  of  the  transition  probability 
matrix . 

Ignall  and  Silver  [1977]  base  their  approach  upon 
results  obtained  by  Buzacott  [1967,  1968]  for  the  limiting 
cases  of  a two-station  production  line  with  no  buffers  and 
with  an  infinite  number  of  buffers.  An  approximation 
formula  for  a two-station  production  line  is  developed  such 
that  the  throughput  rate  is  a linear  combination  of  the 
throughput  rates  of  the  two  limiting  cases  just  mentioned. 

Muth  and  Yeralan  [1981]  discover  a particular  Markov 
model  of  a two-station  production  line  with  an  arbitrary 
number  of  buffers  whose  steady-state  probabilities  possess 
the  scalar-geometric  property.  Using  this  property,  these 
probabilities,  and  therefore  the  throughput  rate,  can  be 
computed  easily.  Yeralan  [1983]  also  shows  that  the  steady- 
state  probabilities  of  a Markov  model  of  any  two-station 
production  line  with  an  arbitrary  number  of  buffers  possess 


the  matrix-geometric  property  provided  the  states  are 
ordered  in  a certain  way. 
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Gershwin  and  Schick  [1983]  make  the  same  assumption 
as  Okamura  and  Yamashima.  Using  a Markov  model  they  calcu- 
late the  steady-state  probabilities  by  assuming  that  these 
probabilities  are  in  a sum  of  products  form.  Their  method 
is  restricted  to  two-  and  three-station  production  lines. 

Shantikumar  and  Tien  [1983]  present  a Markov  model  of 
a two— station  production  line  which  includes  the  possibility 
of  scrapping  an  item  in  a broken  down  station.  Using  a 
concept  called  "level  crossing  analysis,"  a recursive 
solution  procedure  for  the  steady-state  probabilities  is 
given. 

1.2.3.  Models  of  Production  Lines  with  Station  Breakdown 
and  Random  Service  Times 

Buzacott  [1972]  considers  a balanced  two-station  pro- 
duction line.  He  assumes  that  service  times,  repair  times 
and  time  between  breakdowns  are  exponentially  distributed. 
These  assumptions  give  rise  to  a continuous  parameter  Markov 
process,  which  is  solved  for  the  mean  interdeparture  time 
of  items  at  station  2.  Buzacott  compares  exact  results 
obtained  using  this  method  with  results  obtained  using  a 
heuristic  procedure  for  two-station  production  lines,  and 
extrapolates  this  procedure  to  production  lines  with  more 
than  two  stations. 
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Muth  [1979a]  studies  general  properties  of  production 
lines  subject  to  breakdown  without  making  a specific  assump- 
tion on  the  distributions  of  service  times,  repair  times 
or  time  between  breakdowns.  Several  models  for  the  line 
availability  using  different  failure  mechanism  are  given. 
Availability  is  defined  as  the  ratio  of  the  throughput  rate 
of  a production  line  with  station  breakdown  to  the 
throughput  rate  of  the  same  production  line  without  station 
breakdown.  Muth  also  shows  that  the  effects  of  service  time 
variability  and  station  breakdown  can  be  combined  under 
certain  assumptions.  This  means  that  a production  line 
subject  to  breakdown  can  be  represented  by  an  equivalent 
production  line  not  subject  to  breakdown. 

Gershwin  and  Berman  [1981]  present  the  same  Markov 
model  of  a two-station  production  line  as  Buzacott  [1967]  . 
They  calculate  the  steady-state  probabilities  by  assuming 
that  these  probabilities  are  in  a sum  of  products  form. 

They  find  that  these  probabilities  can  be  obtained  by 
solving  a system  of  nonlinear  equations.  The  number  of 
these  equations  is  always  five  irrespective  of  the  number 
of  buffers. 

Berman  [1982]  uses  the  assumption  that  service  on  an 
item  starts  all  over  after  a station  repair.  Through  this 
assumption  he  can  extend  the  result  of  Gershwin  and  Berman 
[1981]  to  the  case  of  Erlang  service  times. 

Altiok  and  Stidham  [1983]  incorporate  the  breakdown- 
repair  phenomenon  into  the  service  time  distribution  as  in 
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Muth  [1979a].  They  assume  that  service  times,  repair  times 
and  time  between  breakdowns  are  exponential.  Using  a Markov 
model,  they  calculate  the  steady-state  probabilities  of  a 
three-station  production  line. 

1.3.  Current  Work 

Even  though  there  has  been  a lot  of  research  conducted 
in  the  production  line  area,  there  are  two  major  features 
that  set  this  dissertation  aside  from  previous  research: 

- The  analysis  is  carried  out  for  an  arbitrary  number 
of  stations  and  for  arbitrary  distribution  functions 
of  service  times,  repair  times  and  time  between 
breakdowns  at  a station. 

- The  numerical  procedures  developed  are  for  a general 
class  of  distribution  functions. 

This  dissertation  is  organized  as  follows.  In  Chap- 
ter 2 we  deal  with  a production  line  model  in  which  station 
breakdown  is  not  considered.  For  a special  class  of  service 
time  distributions  that  we  call  special  phase  type  distribu- 
tions, we  develop  two  numerical  procedures,  one  exact  and 
the  other  approximate,  for  calculating  the  line  throughput 
rate.  We  use  these  numerical  procedures  to  investigate  the 
effect  of  the  number  of  stations  and  of  the  service  time 
variabilities  on  the  throughput  rate  of  a balanced  produc- 
tion line. 

In  Chapter  3 we  incorporate  station  breakdown  and 
repair  in  the  model  of  Chapter  2.  We  derive  a method  to 
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obtain  an  equivalent  production  line  without  breakdown  and 
we  state  conditions  under  which  the  solution  procedures  of 
the  previous  chapter  can  still  be  used. 

In  Chapter  4,  using  the  method  of  Chapter  2,  we  inves- 
tigate the  effect  of  line  unbalance  on  the  throughput  rate. 
We  express  the  throughput  rate  as  a function  of  mean  value 
unbalance  and  variability  unbalance.  We  present  the  results 
in  the  form  of  contours  of  constant  throughput  rate. 

In  Appendix  A we  give  several  properties  of  the  class 
of  special  phase  type  distributions  that  are  important  for 
the  development  of  solution  procedures  of  Chapter  2.  In 
Appendix  B we  specialize  the  method  of  Chapter  2 to  a 
formula  for  the  throughput  rate  of  a three-station  produc- 
tion line.  We  use  this  formula  in  Chapter  4 to  investigate 
the  effect  of  line  unbalance  on  the  throughput  rate. 


CHAPTER  2 


PRODUCTION  LINES  WITHOUT  STATION  BREAKDOWN 
2.1.  Introduction 

In  the  past  production  lines  with  stochastic  servers 
have  typically  been  modeled  as  discrete  state  Markov  pro- 
cesses. In  such  models,  the  service  time  of  the  work 
stations  is  exponentially  distributed.  Extensions  to  Erlang 
service  times  are  feasible  at  the  expense  of  a considerable 
growth  in  the  numbers  of  system  states.  In  Muth  [1977, 

1984]  a holding  time  model  is  presented  that  is  not  based  on 
the  Markov  property  and  is  thus,  in  principle,  not 
restricted  to  any  particular  type  of  service  time  distri- 
bution. Furthermore,  that  model  does  not  suffer  from  the 
explosive  growth  of  the  state  space  as  experienced  by  the 
Markov  model.  The  holding  time  model  leads  to  an  integral 
equation  defining  the  distribution  of  the  interdemand  time, 
and  thus  defining  the  system  throughput  rate.  Closed  form 
solutions  for  the  production  rate  when  there  are  up  to 
three  work  stations,  and  for  special  distributions,  are 
given  in  Muth  [1984] . Numerical  solutions  for  up  to 
three  work  stations  and  more  general  distributions  present 
no  problem.  In  this  chapter  we  extend  the  result  of  Muth 
[1984]  by  developing  a solution  for  the  case  of  more  than 
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three  work  stations.  We  present  two  numerical  procedures, 
one  exact  the  other  approximate,  for  solving  the  integral 
equations  when  the  service  time  distributions  are  of  special 
phase  type.  The  approximate  method  gives  a very  tight  lower 
bound  to  the  actual  production  rate.  It  is  of  interest 
because  it  is  less  complex  than  the  exact  procedure  and 
because  it  can  handle  a larger  class  of  service  time 
distributions . 

The  chapter  is  organized  as  follows.  In  section  2.2  we 
give  the  general  assumptions  pertaining  to  a system  charac- 
terized by  its  states  and  in  section  2.3  we  discuss  the 
corresponding  holding  time  model.  Subsequently  we  establish 
a system  of  integral  equations  for  the  interdemand  time,  and 
we  present  an  exact  numerical  solution  method  for  the  case 
of  exponentially  distributed  station  service  times.  This  is 
followed  by  an  approximate  numerical  method  for  the  same 
case.  Finally,  in  section  2.7  we  generalize  the  numerical 
solution  methods  to  the  case  of  special  phase  type  station 
service  times. 


2.2.  Model  Assumptions 

We  are  treating  a general  model  of  a production  line, 
consisting  of  K work  stations  (servers)  labeled  1,  2,  ..., 
K,  and  arranged  in  series  in  that  order.  There  is  an 
unlimited  supply  of  workpieces,  called  production  items,  at 
station  1 . Each  item  enters  the  line  at  station  1 , passes 
through  all  stations  in  order,  and  leaves  station  K in 
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finished  form.  At  each  station  a service  is  performed.  The 
service  time  of  item  n at  station  j is  denoted  (n) . The 
sequence  (n)  is  a sequence  of  independent  nonnegative 
random  variables,  identically  distributed  as  S ^ . The  random 
variables  S^,  S2,  •••,  are  independent  too;  their  dis- 
tributions are  not  necessarily  identical.  It  is  assumed 
that  a station  cannot  break  down  and  that  each  station  can 
service  only  one  item  at  a time.  Every  station  is  at  any 
time  in  one  of  three  possible  states.  A station  is  busy 
when  it  is  servicing  an  item,  it  is  blocked  when  the  item  on 
which  service  has  been  completed  cannot  move  on  to  the 
subsequent  station  because  that  station  is  busy  or  blocked, 
and  it  is  starved  or  empty  otherwise.  The  starved  state  is 
caused  when  an  item  has  left  the  station,  and  service  on  the 
item  in  the  preceding  station  has  not  yet  been  completed. 

The  durations  of  these  states  associated  with  the  passage  of 
a single  item  are  called  busy  period,  blocking  period  and 
starving  period.  Service  begins  immediately  when  an  item 
moves  into  a station.  Station  1 can  never  be  starved  and 
station  K can  never  be  blocked.  There  may  be  buffer  spaces 
(in-process  storage)  installed  between  stations  to  diminish 
the  occurrence  of  blocking  and  starving.  This  case  does  not 
require  special  treatment  since  a buffer  with  a capacity  of 
L items  is  equivalent  to  L stations  with  zero  service  time 


in  series. 
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With  these  assumptions  the  production  line  is  com- 
pletely specified  given  the  number  K and  the  distribution 
function  F (t)  of  the  station  service  time  S.. 


2.3.  The  Holding  Time  Model 
We  give  a brief  review  of  the  main  features  of  the 
holding  time  model  presented  in  Muth  [1984].  The  model  is 
based  on  a stochastic  process  whose  transitions  occur  when 
an  item  moves  from  one  station  to  the  next.  A move  of  an 
item  from  station  j to  station  j+1  corresponds  to  an  arrival 
at  station  j + 1.  For  je<,  ic  = {1,  2,  ...,  K}  we  define  the 
random  variables 

Bj (n)  = blocking  period  of  item  n at  station  j 
Hj (n)  = period  spent  by  item  n in  station  j,  or 
holding  period  of  item  n at  station  j 
Ij (n)  = starving  period  of  station  j following  the 
departure  of  item  n-1 

Sj (n)  = service  period  of  item  n at  station  j 
and  the  random  vectors 


H(n) 

= [H^(n), 

H2(n),  . 

. Hj^(n)] 

I (n) 

c 

1 — 1 
H 

II 

I2 (n) , . 

• • / ( n ) ] 

S{n) 

= [s^ (n) , 

S2 (n) , . 

• • / ( n ) ] 

From  these  definitions,  we  have 


(n)  = Sj  (n)  + (n) 


(2-3-1) 


20 


The  service  period  process  {S(n);  neN}  is  a superposi- 
tion of  K independent  and  nonidentical  renewal  processes. 

It  induces  the  holding  period  process  {H(n):  neN}  and  the 
starving  period  process  {I(n):  neN}.  The  following 
recursive  relationships  were  derived  in  Muth  [1984]  : 

(n)  = max  [S^  (n)  , Hj_|_^(n-1)  - (n)  ] for  je<  (2-3-2) 

Ij (n)  = max  [Ij_^(n)  + (n)  - Hj(n-l),  0] 

for  jex  (2-3-3) 

To  make  (2-3-2)  and  (2-3-3)  hold  in  all  cases,  (n) , (n) 

and  Sj (n)  are  defined  to  be  identically  zero  for  j = 0 and 
j = K+1.  We  also  have 

Hj^(n)  = Sj^(n)  (2-3-4) 

(n)  = 0 (2-3-5) 

The  system  of  recursive  equations  (2-3-2)  to  (2-3-5)  defines 

how  H(n)  and  I(n)  are  induced  by  S (n) . Of  special  interest 
in  (2-3-2)  and  (2-3-3)  is  the  random  variables 


'Hj (n)  - Ij_^(n+1)  for  je< 


Rj  (n)  = 


0 


otherwise 


(2-3-6) 
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which  represents  the  time  difference  between  the  departure 
of  item  n from  station  j and  the  arrival  of  item  n+1  at 
station  j~l.  This  random  variable  can  be  either  negative  or 
positive.  We  therefore  define 

Rj(n)  = max  [R^  (n)  , 0]  (2-3-6a) 

R~(n)  = min  [R^ (n) , 0]  (2-3-6b) 

The  random  variable  Rj(n)  is  the  residual  holding  period  of 
item  n at  station  j following  the  entrance  of  item  n+1  into 
station  j~l.  The  absolute  value  of  R^ (n)  is  the  time  period 

during  which  both  stations  j and  j-1  are  starved  together. 

Substituting  (2-3-6)  into  (2-3-2)  and  (2-3-3)  we  obtain  for 
all  jeK 

(n)  = max  [S^  (n)  , Rj_|_^(n-1)]  (2-3-7) 

Ij(n)  = max  [Sj_^(n)  - Rj(n-l),  0].  (2-3-8) 

Equations  (2-3-4)  through  (2-3-8)  serve  as  the  basic  equa- 
tions used  to  express  a relationship  between  the  distribu- 
tion functions  of  service  period,  holding  period  and 
starving  period. 

The  relationship  among  the  above  variables  can  be 
nicely  shown  through  activity  networks  as  is  done  in  Muth 


[1984] . The  networks  are  a useful  tool  for  deriving  such 
equations  as  are  used  in  this  section. 
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2.4.  Integral  Equations  for  the  Holding  Time  Model 
If  all  service  times  have  finite  mean  values,  then  it 
is  plausible  that  (n) , (n) , Rt(n)  and  (n)  converge  in 

distribution  to  random  variables  , I ^ , Rt  and  R~  as  n -►  <». 
For  a discussion  see  Muth  [1984].  Thus  from  equations 
(2-3-6)  through  (2-3-8)  the  recursive  equations  for  the 
distribution  function  F of  the  above  quantities  can  be 
obtained. 

If  Hj(n-l)  is  statistically  independent  of  1^  ^ (n)  then 
from  equations  (2-3-6)  and  (2-3-6a)  , for  all  jeic  we  obtain 


3 


/ F (x)dF„  (x+t) 
x=0  ^j-1  «j 


for  t ^ 0 


otherwise  (2-4-1) 


and  from  equation  (2-3-6)  and  (2-3-6b)  , for  all  jei<  we 
obtain 


/ F (x-t)dF„  (X) 
x=0  ^j-1 


for  t < 0 


0 


otherwise  (2-4-2) 
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If  is  statistically  independent  of  (n)  , 

then  from  equations  (2-3-7)  and  (2-3-8),  for  all  je<  we 
obtain 


F (t)  = F (t)F_+  (t) 


(2-4-3) 


0 

1 - / F -(x)dF„  (x+t) 

x=-t  Sj-1 


=< 

3 


/ F_+(x)dF_  (x+t) 

x=0  "j  ®j-l 

for  t ^ 0 


0 otherwise  (2-4-4) 

(see  Appendix  C for  the  derivation  of  these  equations) . 


Equations  (2-4-1)  through  (2-4-4)  are  the  system  of 

equations  defining  F (t) , F^  (t) , F_+(t)  and  F„-(t)  for  all 

H.  I.  R.  R, 

jeic.  For  j - 0 and  j = K+1  they  are  equal  to  u(t)  which  is 
the  unit  step  function  with  jump  point  at  t = 0. 

It  follows  from  the  assumed  statistical  independence  of 
Sj  (n)  that  (n)  and  (n)  are  also  statis- 

tically independent,  which  can  be  seen  in  the  networlc 
representation  given  in  Muth  [1984]  . The  question  of 
whether  or  not  (n-1)  and  I^_^(n)  are  statistically  inde- 
pendent when  K > 3 is  a difficult  one.  We  have  not  been 
able  to  resolve  this  question  from  an  analysis  of  the 
complex  relationships.  Therefore  we  have  conducted  exten- 
sive simulations  of  four-  and  five-station  balanced  produc- 
tion lines  with  exponential  service  times.  As  a result  of 
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100  simulation  runs,  the  95%  confidence  interval  for  the 
sample  correlation  coefficient  of  Hj(n-l)  and  ^(n)  is 
[0.003,  0.006].  This  result  indicates  a very  slight  posi- 
tive correlation.  For  practical  purposes,  however,  we 
may  conclude  that  (n-1)  and  Ij_^(n)  are  statistically 
independent . 

Even  when  the  index  j represents  a station  and  the 
index  j-1  represents  a buffer,  or  vice  versa,  the  sample 
correlation  coefficient  of  (n-1)  and  Ij_^(n)  is  about  the 
same  as  above.  However,  when  the  indices  j-1  and  j are  both 
for  buffers,  the  sample  correlation  coefficient  of  (n-1) 
and  Ij_^(n)  is  much  higher.  For  two-station  balanced 
production  lines  with  a buffer  of  capacity  two  or  three  and 
exponential  service  times,  the  95%  confidence  interval  for 
the  sample  correlation  coefficient  of  Hj(n-l)  and  ^ (n)  is 
[0.03,  0.09]. 

Substitution  of  (2-4-3)  into  (2-4-1)  and  (2-4-2)  gives, 
for  all  jeK, 


FRt(t)  = 
3 


1 - f F (x)d(F  (x+t)F  + (x+t)  ) for  t ^ 0 

x=0  -^j-1  3 + 1 


otherwise  (2-4-5) 


FR-(t)  = 
3 


1 - / F (x-t)d(F  (x)F  + (x)  ) for  t < 0 

x=0  -^j-1  ^j  + 1 

otherwise  (2-4-6) 
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Given  the  set  of  F (t)  and  F (t)  for  all  jeK,  the  set 

Ij  Sj 

of  Fj^+(t)  can  be  generated  using  equation  (2-4-5)  and 

j 

backward  substitution  starting  with  F (t)  = u(t). 

^K+1 

Likewise  given  the  set  of  F„+(t)  and  F_  (t)  for  all  je<,  the 

R.  S, 

set  of  F^-(t)  and  F (t)  can  be  generated  alternately  using 

j j 

equation  (2-4-4)  and  (2-4-6)  and  forward  substitution 

starting  with  F (t)  = u(t)  or  F_- (t)  = u(t)  (since  I,  s 0, 

1 ^2  ^ 

R2  = H2  which  is  always  positive) . Thus  we  have  three 

functional  relationships  between  the  set  of  F_+(t),  the  set 

R « 

3 

F -(t)  and  the  set  of  F (t) . In  combination  these  three 

i . . . 3 

relationships  establish  a system  of  simultaneous  equations 

defining  F„+(t),  F_-(t)  and  F_  (t) . 

R.  R.  I, 


2.5.  The  System  Throughput  Rate 
The  throughput  rate  r of  the  production  line  is  defined 
as  the  number  of  items  entering  the  first  station  per  unit 
time,  in  the  long  run.  Let  N(t)  be  the  number  of  items 
entering  the  first  station  in  the  time  interval  (0,t).  Then 


r 


lim 

t->-oo 


E[N(t)  ] 
t 


(2-5-1) 


This  is  equivalent  to 


1 

E[H^] 


r 


(2-5-2) 
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where 

00 

E[HJ  = / (1-F„  (t))dt  (2-5-3) 

The  random  variable  is  the  interdemand  time  of  station  1. 
Since  the  first  station  is  never  starved,  the  interarrival 
time  of  items  to  the  production  line  is  distributed  as  . 

By  the  assumption  that  (n)  converges  in  distribution  to 
H j , has  a finite  mean.  The  equivalence  between  (2-5-1) 

and  (2-5-2)  thus  follows  from  the  law  of  large  number  (see 
Cox  [1962]).  From  (2-4-3)  it  follows  that 

F (t)  = F (t)F  +(t)  (2-5-4) 

1 ^1  ^2 


Since  the  first  station  is  never  starved  it  is  seen  from 
(2-3-6)  that  R2  = H2;  hence  using  (2-4-3)  again  we  have 

F (t)  = F (t)F  (t)F_+(t)  (2-5-5) 

1 ^1  ^2  ^3 

To  compute  F„  (t)  we  determine  F_+(t)  from  the  simultaneous 

Hi  R3 

equations  established  by  (2-4-4)  through  (2-4-6) . A solu- 
tion procedure  for  the  case  where  S^  for  all  je<  are  expon- 
entially distributed  is  presented  in  the  following  section. 

2.6.  Solution  for  Exponential  Service  Times 
We  are  treating  the  case  where  all  stations  have 


exponential  service  times,  i.e. 
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-a  . t 


1-e 


for  t & 0 


F„  (t)  = 


(2-6-1) 


otherwise 


For  this  choice  of  service  time  distribution  F + (t)  takes 

"j+i 

the  following  form  (see  appendix  A for  a formal  proof) : 


F„+  (t) 

"j+i 


m . 

3 

Z a.  . 
i=l 


for  t ^ 0 


otherwise 


(2-6-2) 


where  s . . >0  and  OS  Z a . . S 1 . 

i=l 

Equation  (2-6-2)  can  be  written  as 

m . 

3 -s..t 

F + (t)  = Z c.  .e  for  t > 0 (2-6-3) 

^j+1  i=0 


where  c^^  = 1,  s^^  = 0 and  c^^  = -a^^  for  i ^ 0. 
Substitution  of  (2-6-3)  into  (2-4-5)  gives 

m . 

3 

F +(t)  = 1 + Z c.  .G(s.  . , j ,t) 
i=0 


(2-6-4) 
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where 


-s . . (x+t) 


(s^.,j,t)  = - / Fj  (x)d{Fg  (x+t)e  ) 


x=0  "j-1 


-s . . t 
ID 


= A(j-l,s.  .)e  - A{j-l,s.  .+a.)e 


- ( s . . +a . ) t 
ID  D 


ID 


ID  D 


(2-6-5) 


and 


A(  j ,s) 


00 

( / F (x) se  ®^dx 

jx=0 

0 


for  s 7^  0 


otherwise 


(2-6-6) 


Equations  (2-6-2)  through  (2-6-6)  are  the  recursive  equa- 
tions that  can  be  used  to  generate  F_+(t)  from  F_+  (t) . 

j j + 1 

However,  one  problem  with  these  equations  is  that  the 

parameters  a^^^  and  s^^  are  still  unknown.  Therefore  we  will 

bring  the  equations  into  a form  which  does  not  contain  a. . 

ID 

and  Sj^ j . Toward  this  end  we  define  a linear  operator  * as 
follows : 


Fg  (t)  * e"®^  = A(j-l,s)e"®^ 


- ( s+a  . ) t 

A(j-l,s+Oj)e 


(2-6-7) 
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By  comparing  the  right  hand  sides  of  (2-6-5)  and  (2-6-7) , 
and  noting  the  form  of  F + (t)  it  is  easy  to  see  that 

equation  (2-6-4)  can  be  rewritten  as 


F +(t)  = 1 + F (t)  * F_+  (t) 

^j+1 


(2-6-8) 


Equation  (2-6-8)  represents  a method  that  recursively 

generates  F +(t)  from  F_+  (t)  starting  with 

J ^j+1 

(t)  = e^^  E 1.  For  j=K  and  j=K-l,  the  explicit  form  of 

K+1 

F +(t)  is  shown  below. 


F +(t) 


= 1 + F„  (t) 


Ot 


= 1 + A(K-l,0)e°^  - A(K-l,a^)e 

K 


~“k^ 


= 1 - A(K-l,ctj^)e 


-CRt 


F + (t)  = 1 + F„  (t)  * (l-A(K-l,a^)e  '"^^) 

^-1  ^K-1  ^ 


= 1 - A (K-2  , J^)  e - A(K-l,Oj^)  X 


(A(K-2 


- A(K-2,aj^+a^^_^) 
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To  make  the  procedure  useful  for  numerical  calculation, 
define  a function  Gj^(j,s)  as  follows: 


Gq  ( j ,s) 

-St 

= e 

G^ (j ,s) 

= F (t)  * e"®^ 

G2 ( j f s) 

= F (t)  * F (t)  * e"®^  = F„  (t)  * G^(j+l,s) 

Sj  S.  1 

Gj^(j  ,s) 

= F (t)  * F„  (t)  * ...  * F„  (t)  * e"®^ 

^j+1  ^j+k-1 

(2-6-9) 

Gj^(jfS)  is  a superposition  of  exponential  functions.  The 
number  of  exponential  terms  in  Gj^(j,s)  is  2^.  From 
equations  (2-6-7)  and  (2-6-9)  we  can  see  that  Gj^(j,s) 
satisfies  the  following  recursion: 


Gq ( j f s) 

-st 

= e 

Gj^(  j ,s) 

= A(  j+k-2,s)Gj^_^  ( j ,s) 

- A(j+k-2,s+a.^j^_^)  G3^_^(j,s+a.^j^_^) 

for  k = 1 , 2 , ...  (2-6-10) 
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Using  the  factor  defined  below  and  applying  equa- 
tions (2-6-7)  through  (2-6-9),  the  evolution  of  F +(t)  is 

R • 

3 

shown  for  several  values  of  j as  follows: 


j 

j 


and  so  on. 

Thus,  noting  that  GQ(j,s)  = e and  e^^  = 1 , we  have  that 


K 

F +(t)  = 1 + E c.G.  .(j,cx.)  (2-6-11) 

Rj  1 1 


where,  from  equation  (2-6-7)  with  s = 0, 


= -A(i-l,a^) 


c . 
1 
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Hence 


K 

Fj^+(t)  = 1 - Z A(i-l,a.)G.  .(j,a.)  (2-6-12) 

j i=j  -> 

and  we  have  the  means  to  compute  Gj^(j,s)  readily  at  hand 
given  by  equation  (2-6-10) . As  an  illustration,  we  give  the 
expressions  of  F^^+(t)  for  j=K  and  j=K-l  below. 

j 


F +(t) 


1 - A{K-1,cl^)Gq(K,cl^) 


= 1 - A(K-l,aj^)e 


Fr+  (t)  = 1 - A(K-2,Or_^)Gq (K-1,Or_^)  - A(K-l,aR)G^(K-l,a  ) 

K—  1 


= 1 


A(K-2,aR_^)e 


““k-1^ 


A(K-l,aR)  X 


(A(K-2,Or) 


- A(K-2,OR+aR_^) 


As  a next  step,  it  is  necessary  to  establish  a recursive 
expression  for  (t)  which  corresponds  to  the  expressions 

j 

for  F + (t)  given  in  (2-4-3) . From  equation  (2-4-3)  and 

j + 1 

(2-6-12)  we  have 
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3 


3 


K 

E A(i-l,a. )G' . . . (j+l,a. ) 

i=j+l  ^ i-j-i  1 


(2-6-13) 


where 


Gj^(j+l,s)  = Gj^(j  + l,s)Fg  (t)  . 


(2-6-13a) 


By  comparing  (2-6-13a)  with  the  recursion  (2-6-10)  we  see 
that  Gj^(j  + l,s)  satisfies  that  same  recursion  (2-6-10) 
given  for  Gj^(j+l,s)  except  that  it  starts  with 
G'  (j  + l,s)  = e"®^Fg  (t)  . 

j 

We  define  the  parameter  as 


/ F +(x)a.  ^e““j-l^dx 

0 


for  j=2 , 3 , . . . ,K 


otherwise  (2-6-14) 


Substituting  (2-6-12)  into  (2-6-14)  gives 


B . 


3 


1 


K 

E A(i-l,a^)E^_ . (j,a^)  forj=2,3 


,K 

(2-6-15) 


where 
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EK(jfS)  = / ( j , s) aj_^e““j-l^dx; 


E^(j,s)  satisfies  the  recursion  (2-6-10)  given  for  Gj^(j,s) 
except  that  it  starts  with 


Eo(j/S)  = Gq ( j ,s) aj_^e"“j-l^dx  = aj_^/ (s+aj_^)  . 


At  this  point  we  have  completed  the  derivation  of 
recursive  equations  for  Fj^+(t)  and  F (t)  . The  next  step  is 

j j 

to  derive  similar  equations  for  F (t)  and  F -(t).  However 

I.  R. 

F^  (t)  and  F -(t)  are  not  of  special  phase  type  (see 

j j 

appendix  A) . It  turns  out  that  the  form  of  F^  (t)  and 

Fj^-(t)  depends  on  the  particular  condition  of  the  production 

line.  An  exact  expression  for  F (t)  and  F -(t)  for  the 

I.  R. 

balanced  case  is  obtained  next.  The  unbalanced  case  has 
little  practical  interest. 


2.6.1.  The  Exact  Method;  Balanced  Line 

The  derivation  so  far  allowed  for  the  station  service 

times  to  be  nonidentical.  From  now  on,  in  order  to  be  able 

to  keep  the  derivation  of  F (t)  and  F -(t)  simple,  we 

Ij  R, 

require  identical  service  times,  and  we  therefore  put 
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We  alternately  use  equations  (2-4-4)  and  (2-4-2)  to 


derive  F (t)  and  F - (t) , starting  with  j=2  and  moving 

j j + 1 

forward  toward  j=K-l.  From  = 0,  it  follows  that  R2  = 
so  that  R2  is  always  positive.  Hence  from  equations  (2-4-4) 
and  (2-6-14) 

Fj  (t)  = 1 - £26"“^ 


where  C is  some  constant.  From  equation  (2-4-4)  again 


2 


and  from  (2-4-2) 


F -(t) 
^3 


B2Ce 


at 


(t)  = 1 
3 


x=-t 


-at 


-at 


B2Cate 


-at 


Continuing  in  this  fashion,  we  obtain 


j 


(2-6-16) 
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Substitution  of  (2-6-16)  into  (2-4-2)  gives 


j-1 


oo  i_l 

a 


F - (t)  = Z D ( / - , (x-t)^"^e"“''dF  (x))e“^ 

:+l  i=l  x=0  • “j+l 


j-1  i-1  i-1  / “>  . . . , ^ 

^ ^ii  ^ ^ in — M-k-1  ) I ^ s dF  (x))t  e 

i=l  k=0  ^ • x=0 


j + 1 


j-1  i-1  / „vk  “ I „.i-k-l  , 

T*  r\  T'  (""Ot)  / f (otX)  / \ \ I k Ott 

°ij  v!„  *„£„  ' d-k-l)l  ® ® 


i=l  k=0  x=0 


j-1  i-1  , ..k 

E D.  . E C . . , € 

i=l  k=0 


at 


(2-6-17) 


where 


OO  i “ 1 

„ _ f (ax)-'  -aXj-  , V ^ 

i.j  ■ 0 ■<]-!)  ‘ ® “i+l  3 < 1 


(2-6-18) 


and  where  F^^  (x)  is  given  by  (2-6-13)  . Substituting 

j 

(2-6-13)  into  (2-6-18)  we  obtain 


K 


C.  . = 
13 


— -j-  + E A(k-l,a)E'  . _(i+2,j,a)  for  j < i 

2^  k=i+2  K-i-z 


otherwise  (2-6-19) 


where 
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(ax) ^-ax 


Ej^(i,j,s)  satisfies  the  same  recursion  as  Gj^(i,s)  which  in 
turn  follows  the  same  recursion  as  Gj^(i,s).  The  starting 
value  is 


EJ(i,j,s)  = r e-“='aG'(i,s)  = 


(s+2a) ^ (s+a) ^ 


In  order  to  find  F (t)  we  evaluate  the  following 

j+1 

integral: 


0 „ j-1  i-1  / 0 , 

/ F„-  (x)ae"“  dx  = E D,  ^ E C.  . / ax^dx 


-t  ^j+1 


i=l  k=0 


k!  j/i-k  , 

' x=-t 


j-1  i-1 

= E D.  . E C . . 


(at) 


k+1 


i=l  k=o 


Thus  from  (2-4-4)  and  (2-6-14)  we  obtain 


. j-1  i-1  , ,>k+l 

F (t)  = 1 - B.  ^e"“  - E D.  . E C.  . , - . , e““^ 

Ij  + 1 3 + 1 13  3,1-k  (k+D! 

(2-6-20) 


-ED 

i=l 


i, j+1  (i-1) ! 


(at)^  ^ »~at 
e 


1 


(2-6-21) 
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i-1 


Comparing  the  coefficient  of  , e““^  in  both  (2-6-20) 

and  (2-6-21) , we  have  the  following  recursive  equations  for 


°i,j+l'  i 


°i, j+1 


°k,j  ^j,k-i+2 


for  i = 1 

(2-6-22) 

for  i = 2 , 3 , . . . , j 


Finally,  substitution  of  (2-6-16)  into  (2-6-6)  gives,  for 
s 0,  ' 


j-1 

1 - Z D.  . 
i=l 


A(j,s)  = 


i-1 
a ^ 

(s+a) ■ 


for  j = 2 , . . . , K 


otherwise 


(2-6-23) 


Equations  (2-6-15),  (2-6-19),  (2-6-22)  and  (2-6-23) 
constitute  the  system  of  equations  that  needs  to  be  solved 
in  order  to  find  the  value  of  i < j,  j = 2,3,...,K-1 

is  not  needed)  . Once  the  have  been  determined,  the 

system  throughput  rate  can  be  calculated  by  substituting 
equation  (2-6-12)  for  j=3  into  (2-5-5)  to  obtain  the  mean 
interdemand  time;  that  is 


K 

E[H^ ] = y + Z A(i-l,a. )E"  - (a. ) 

i=3  ^ ^ 


(2-6-24) 
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where 


00 

p = / (1-F  (t)F  (t))dt  (2-6-25) 

0 ^2 


E^(s) 


/ Gj^(3,s)Fg  (t)Fg  (t)dt; 


(2-6-26) 


E^(s)  satisfies  the  recursion  (2-6-10)  given  for  Gj^(3,s) 
except  that  it  starts  with 


00 

E"(s)  = / e"®^F  (t)F„  (t)dt. 

0 ^2 


The  number  of  equations  that  need  to  be  solved  is  in 
the  order  of  K . The  equations  are  nonlinear  and  are  in  the 
form  of  a fixed  point  problem,  i.e.  the  problem  of  finding  a 
vector  X which  satisfies  x = f (x)  where  f is  a vector 
valued  function.  We  know  that  because  of  statistical  equi— 
and  are  unique  and  all  of  their  components 
must  have  values  between  0 and  1.  This  information  helps  in 
choosing  a good  initial  value  for  the  solution  procedure.  A 
simple  iterative  scheme  can  be  devised  as  follows: 

Step  1:  choose  and  let  k=0 . 

Step  2:  calculate  using  (2-6-22) 

calculate  using  (2-6-15)  and  using 


(2-6-19) 
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Step  3:  if  - B ^ | < e and  " ^ij  I ^ ^ 

all  i and  j,  stop. 

Otherwise  let  k = k+1  and  return  to  step  2. 
This  method  has  been  found  to  work  well. 

It  is  instructive  to  compare  the  method  presented  with 
that  of  solving  a Markovian  state  model.  Both  involve  the 
solution  of  a system  of  equations.  The  Markov  model  of  a K 
station  line  has  a total  number  of  states  which  grows 
asymptotically  as  (2 . 62) ^/2 . 24 . Thus  the  case  K=6,  which  is 
the  largest  value  of  K for  which  a solution  is  given  in  the 
literature,  results  in  144  states  (it  grows  to  6765  for 
K=10) . This  means  that  a system  of  144  simultaneous 
equations  must  be  solved  to  obtain  the  steady  state 
probabilities.  In  the  method  of  this  paper  there  are 
(K-2)(K-l)/2  simultaneous  equations  (i.e.  the  number  of 
variables  B^  and  that  are  needed  to  calculate  the  system 
throughput  rate)  to  be  solved.  The  convergence  of  the 
solution  procedure  given  above  has  been  found  to  be  very 
fast.  For  the  case  K=10,  it  took  18  iterations  with 
30  seconds  of  run  time  on  a VAX  11/750  under  the  VAX/VMS 
operating  system.  Throughput  rates  computed  for  values  of  K 
up  to  10  are  shown  in  Figure  2.1  where  r^  = E[S^]/E[H^], 
called  the  utilization  or  the  normalized  throughput  rate,  is 
graphed  as  a function  of  K.  Comparing  the  results  of  this 
procedure  and  the  result  of  Hillier  and  Boling  [1967]  shows 
that  the  difference  in  the  case  K=6  is  about  0.1%.  This 
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small  difference  could  be  due  to  the  numerical  accuracy  of 
either  method. 

2.6.2.  An  Approximation  Method 

Since  represents  the  period  of  time  where  station  j 


probability  of  being  zero,  especially  in  a balanced  line. 
Therefore,  as  an  approximation,  we  may  assume  that  R~  is 
zero  with  probability  one.  In  terms  of  the  defining 
equations  this  assumption  affects  only  the  expression  for 
(n) ; i.e.  equation  (2-3-8)  becomes 

Ij  (n)  = max  [Sj_^  (n)  - Rj(n-1),0]  for  jeic  (2-6-27) 


Substitute  (2-6-28)  into  (2-6-6)  to  obtain,  for  s # 0, 


Hence 


00 


F£  (t) 


1 - / F +(x)dF  (x+t) 

x=0  ®j-l 


for  t > 0 (2-6-28) 


1 - BjS/ (s+ttj_^)  for  j=2,3 


K 


A(  j ,s) 


(2-6-29) 


1 


otherwise 
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The  solution  of  equations  (2-6-29)  and  (2-6-15)  yields  the 
values  of  B ^ , j=2,...,K-l.  Once  the  are  found  the  line 
throughput  rate  can  be  calculated  using  equation  (2-6-24) 
where  A(j,s)  is  defined  as  in  (2-6-29).  In  this 
approximation  only  (K-2)  equations  in  the  form  of  a fixed 
point  problem  need  to  be  solved.  Furthermore  the  equations 
themselves  are  much  simpler.  This  reduction  of  complexity 
results  in  a considerable  reduction  in  computational  effort. 
For  the  case  K=10,  it  took  8 iterations  with  5 seconds  of 
run  time.  Computed  approximations  are  included  in  Fig- 
ure 2.1. 


An  advantage  of  this  approximation  method  is  that  it  is 
applicable  to  a larger  class  of  service  time  distributions. 
Moreover,  the  same  solution  procedure  applies  for  identical 
and  nonidentical  server. 

The  approximate  throughput  rate  is  less  than  the  actual 
throughput  rate  since  is  never  larger  than  and  there- 
fore the  corresponding  is  never  smaller  than  when  it  is 
calculated  using  Ij;  see  equation  (2-3-2) . Therefore,  by 
using  this  approximation,  E[H^]  is  larger  than  its  actual 
value.  It  should  be  noted  that  any  error  introduced  into 
the  distribution  F£  (t)  using  equation  (2-6-28)  will  subse- 

j 


quently  be  diminished  considerably  when  the  corresponding 
F + (t)  is  used  in  equation  (2-4-6) . To  explain  the  dimin- 

j + 1 

ishing  effect  consider  equation  (2-3-7) . The  error  present 
in  the  value  of  any  particular  realization  of  1^ (n)  will 
become  an  error  in  the  corresponding  realization  of  (n) 
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only  if  Rj_^^(n-1)  = (n)  > (n)  , since  other- 

wise Hj  (n)  = Sj (n) . This  condition  will  exist  only  for  a 
small  fraction  of  all  realizations. 


2.7.  Solution  for  Special 
Phase  Type  Service  Times 

Let  the  service  time  distribution  of  the  jth  station  be 
of  special  phase  type  discussed  in  Appendix  A;  i.e. 


Fs  (t)  = 1 - 

1 


m . 

i=l 


-a  . . t 

ifD 


for  t S 0 


(2-7-1) 


Equation  (2-6-8)  still  holds  here,  but  with  the  operation  * 
defined  as  follows: 


Fg  (t)  * e“®^ 

j 


A ( j-1 ,s) e 


_ , m . 

® - Z^A ( j-1 ,a . .+s) e . 

i=l  -‘•'J 


-(s+a.  .)t 

•^  / J 


(2-7-2) 


where  A(j,s)  is  defined  as  in  equation  (2-6-6).  From  the 
definition  of  Gj^(j,s)  as  given  in  equation  (2-6-9),  the 
expression  of  F^+(t)  must  take  the  following  form: 


FR^(t)  = 


K 

- Z 

k 


m. 


for  t 


1 


(2-7-3) 
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and  from  (2-7-2)  Gj^(j,s)  satisfies  the  following  recursive 
relation: 

GQ(j,s)  = e"®^ 

Gk(j/S)  = A(  j+k-2,s)Gk_3^  ( j ,s) 

D+k-l 

i=l  j+k-l^S-1  j+k-1^ 

(2-7-4) 


Similarly,  from  equation  (2-4-3)  and  (2-7-3)  we  have 


F (t)  = F„  (t) 
H.  S. 


K m, 

Z 6.  ,A(k-l 

k=j  i=l 


“i,k^^'k-j-l 

(2-7-5) 


where  G> _ ( jtl , a . = ^k- j-1 < ' “i ,k> ^ 


Therefore  Gk(j+l,s)  satisfies  the  recursion  given  in  (2-7-4) 
except  that  it  starts  with  G^(j+l,s)  = GQ(j+l,s)Fg  (t) . 

The  mean  interdemand  time  is  obtained  by  substituting 
(2-7-3)  for  j=3  into  (2-5-5) 


K 


m. 


E[H^]  = 


y + 


^3  i=l 


(2-7-6) 
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where  p and  E^(s)  are  defined  as  in  (2-6-25)  and  (2-6-26), 
respectively.  This  quantity  can  be  calculated  only  after 
A(k-l,s)  is  found.  To  do  this  we  need  to  define 


“ -a  X 

f F +(x)a  . ^e  dx  for  j=2 

0 j 


. ,K 

(2-7-7) 


Substitute  (2-7-3)  into  (2-7-7)  to  obtain 


K m. 


1 - Z 
k 


=3  i=l  ^'“i,k^\-j  '“i,k^ 

for  j=2 , 3 . . . ,K 


^ 


otherwise 


(2-7-8) 


“ -a  . ^ X 

where  E,  (£,j,s)  = / G,  (j,s)a.  . ,e  dx  ; 


Ej^(^,j/S)  satisfies  the  recursion  given  in  (2-7-4)  except 
that  it  starts  with  E (£,j,s)  = a . ./(a  . ^+s). 


2.7.1.  The  Approximation  Method 

Substitution  of  (2-7-1)  into  (2-6-28)  for  j = 2,3, ,K 

gives 


"'j-1 


F (t)  = 1 - r B.  .8.  . ,e 

Ij  i=l  3/1  1.3-1 


-a . ■ 1 t 

1/3-1 


for  t ^ 0 (2-7-9) 
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Substitution  of  (2-7-9)  into  (2-6-6)  gives,  for  s # 0, 


A(j,s)  = 


^ ■ ifi  = 3=2,3 K 

Otherwise  (2-7-10) 


The  system  of  equations  defined  by  (2-7-8)  and  (2-7-10) 
must  be  solved  to  determine  the  value  of  B . . ; j=2,...,K-l 
and  i=l , 2 , . . . ,mj . The  system  throughput  rate  is  then 
obtained  from  equation  (2-7-6)  where  A(j,s)  is  defined  as  in 
(2-7-10) . The  number  of  equations  that  need  to  be  solved 
K-1 

here  is  I m. , which  means  it  grows  linearly  with  the  number 
i=2  ^ 

of  stages  in  the  service  time  distributions  used.  As  a 
comparison  in  the  Markov  model  this  number  grows  polynomially 
of  degree  K with  the  number  of  stages. 

Using  this  method  we  calculated  the  normalized 
throughput  rate  of  a balanced  line  for  the  case  of  m-stage 
general  Erlang  service  times  for  up  to  K=8  and  m=4.  The 
results  are  shown  in  figure  2.2,  where  r^  is  graphed  as  a 
function  of  rig  (squared  coefficient  of  variation  of  service 
time  S)  with  K as  a parameter.  Since  the  throughput  rate  of 
a balanced  production  line  is  predominately  determined  by  the 
first  two  moments  of  the  service  time  distribution  (Muth 
[1977]),  this  graph  can  be  considered,  for  all  practical 
purposes,  to  be  a distribution-free  property  of  a serial 
production  line. 
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Figure  2.2.  The  throughput  rate  of  K stations  production 

2 

line  as  a function  of  rig:  Erlang  Case. 
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2.7.2.  The  Exact  Method 

Similar  to  what  we  did  for  the  exponential  case,  if  the 
line  is  balanced,  i.e.  if  = a^,  and  = m for 

all  jsK,  we  have 

m 

A(2,s)  = 1 - I B . 3 . s/ (s+a . ) 

i=l  ^ ^ 


m 

A(3,s)  = 1 - Z B .B.s/(s+a.)  - 

i^l  J,i  1 1 


m m 

( Z B .S.C.s  (s+a.))  ( Z 3.a./(s+aJ). 

i=l  1 1 1 3 3 


where 


00  -a . X 

C.  = f e ^ dF„  (X) . 

0 “3 

The  expression  for  A(i,s)  when  i > 3 grows  very  quickly 
in  complexity.  We  are  still  able  to  derive  the  expression 
for  A(4,s),  beyond  that  things  appear  hopelessly  complicated. 
Recall  from  equation  (2-7-6)  that  A(i,s),  i = 2,  . . .,  K-1 
are  required  to  compute  the  throughput  rate  of  a K-station 
production  line.  For  this  reason  the  computer  program  was 
written  to  handle  cases  with  K ^ 4.  The  numerical  results  of 
K = 4 shows  that  for  the  case  rig  = 1//2  the  difference 
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between  the  throughput  rate  values  calculated  using  the  exact 
method  and  the  approximation  method  is  a half  of  the 
difference  for  the  case  rig  = 1.  Since  the  approximation 
method  is  good  enough  for  the  case  rio  = 1/  we  may  say  that 
for  rig  < 1 it  is  sufficient  to  use  the  approximation  method. 

2.8.  An  Empirical  Formula  for  the  Throughput  Rate 
of  a Balanced  Production  Line 

The  throughput  rate  of  a paced  production  line  is  given 
by 


E[max(S^,S2,  . . . ,Sj^)  ] (2-8-1) 

Muth  [1973]  shows  that  (2-8-1)  is  a lower  bound  for  the 
throughput  rate  of  the  unpaced  model  of  the  same  production 
line . 

An  explicit  expression  for  the  lower  bound  (2-8-1)  in 
the  case  of  identical  and  uniformly  distributed  service  times 
with  mean  value  equal  to  one  is  derived  in  Muth  [1973].  The 
expression  is 


r 


L 


1_ 

1 + /3 


K-1 
'^S  K+1 


(2-8-2) 


The  uniform  distribution  was  chosen  because  it  covers  a wide 
range  of  rig;  that  is  0 S rig  ^ 1//3.  However,  formula  (2-8-2) 
is  intended  to  be  used  for  all  values  of  n„. 
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The  relative  difference  between  the  lower  bound  (2-8-2) 
and  the  actual  throughput  rate  is  a monotonically  increasing 
function  of  rig.  This  means  that  the  difference  is  bigger  for 
a bigger  value  of  rig.  We  compare  the  values  of  throughput 
rate  calculated  using  (2-8-2)  with  the  values  of  throughput 
rate  calculated  in  previous  sections  using  general  Erlang 
service  time.  For  K = 3 with  rig  = 1//3,  the  difference  is 
about  2.2%  and  with  rig  = 1 / the  difference  is  about  5%. 

Considering  the  form  of  (2-8-2)  we  assume  that  an 
improved  approximation  for  the  actual  normalized  throughput 
rate  should  be  of  the  form 


where  f(rigfK)  is  a function  with  the  following  properties: 

- f(0,K)  = 0 

- f(ng,l)  = 0 

- nondecreasing  with  respect  to  n and  K. 

We  assume  further  that  f(rig,K)  in  (2-8-3)  has  the  following 
form: 


The  form  of  (2-8-3)  together  with  (2-8-4)  is  very  similar  to 
the  form  (2-8-2)  except  that  there  is  an  additional  term  bn 
in  the  denominator  of  f(no,K). 


r. 


1 


N 1 + f(ng,K) 


(2-8-3) 


K+l+bng 


(2-8-4) 


This  additional  term  serves 
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as  a correction  factor  that  will  compensate  the  difference 
between  the  lower  bound  and  the  actual  value  which  is  bigger 
for  bigger  rig. 

The  task  is  now  to  determine  the  values  of  a and  b in 

(2-8-4)  from  the  numerical  values  of  r„  obtained  by  the 

method  of  previous  sections.  From  each  value  of  r.  we 

N 

computed  f(rig,K)  using  equation  (2-8-3)  and  fitted  these 
values  of  f(rigfK)  into  equation  (2-8-4)  to  get  the  values  of 
a and  b that  minimize  the  error.  We  obtained 
a = 1.67 
b = 0.31. 

Hence 


r 


N 


1 + 1.67 


K-1 

K+1  + . 3 1 rig 


In  this  curve  fitting  we  used  K=3,  4,  . . .,  10  and 

Tig  = 1//3,  1//2,  1.  The  maximum  error  obtained  in  these 

ranges  of  values  of  K and  ri_  is  about  1%.  For  values  of 

b S 

that  are  less  than  1//3,  we  expect  to  get  a smaller  error 
since  this  is  the  case  for  K = 3. 

The  absolute  throughput  rate  r is  obtained  by  dividing 
r^^  with  the  mean  service  times;  that  is 


r 


'N 


/y 


where  y is  the  mean  service  time. 


CHAPTER  3 


PRODUCTION  LINES  WITH  STATION  BREAKDOWN 
3.1.  Introduction 

In  this  chapter  we  treat  the  case  where  stations  are 
subject  to  breakdown  and  repair.  For  this  case,  the  fol- 
lowing assumptions  are  made.  The  number  of  breakdowns  of  a 
station  during  a fixed  time  interval  is  a random  variable. 

A station  can  break  down  only  while  it  is  actively  servicing 
an  item  (i.e.  in  the  busy  state).  When  a station  breaks 
down,  the  service  is  interrupted  and  a repair  is  begun.  The 
item  in  progress  remains  in  the  station.  Service  is  resumed 
on  the  item  in  progress  after  the  repair  has  been  completed. 
Breakdown  of  a station  does  not  cause  a shutdown  of  the 
entire  line,  but  may  affect  adjacent  stations  through 
starving  and  blocking.  Under  the  above  assumptions,  the 
period  of  time  during  which  a station  is  broken  down  will 
always  be  nested  in  a service  period  for  a single  item. 

The  key  idea  of  the  following  development  is  that  the 
total  time  elapsed  from  the  beginning  to  the  completion  of  a 
service,  including  any  repair  periods  that  may  intervene,  is 
treated  as  a single  random  variable.  This  period  (denoted 
by  V)  will  be  called  the  virtual  service  period.  Any 
analysis  method  applicable  to  a model  that  does  not  consider 
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breakdown  and  uses  as  input  the  statistics  of  the  service 

times  may  therefore  be  applied  to  the  model  including 

breakdown,  by  using  the  statistics  of  the  virtual  service 

period  instead  of  the  statistics  of  the  service  times.  More 

specifically,  if  the  throughput  rate  r^  of  the  production 

line  model  without  breakdown  is  some  function  f of  the 

distributions  F (t) , F (t) , . . . , F„  (t)  of  the  station 
^1  ^2 

service  times,  namely 

To  = f [Fg  (t) , Fg  (t) , . . . , Fg  (t) ] 

12  K 

then  the  throughput  rate  r^  of  the  same  model  including 
breakdown  is 

r = f [F  (t)  , F (t)  , . . .,  F..  (t)] 

1 2 

where  F^  (t) ,i=l,  2,  . . .,Kis  the  distribution 
i 

function  of  the  virtual  service  period  of  station  i. 

The  chapter  is  organized  as  follows.  In  section  3 . 2 we 
derive  a general  formula  for  the  transform  of  F^(t).  In 
section  3.3  we  specialize  this  formula  to  several  models  of 
the  breakdown  counting  process.  In  section  3.4  we  consider 
the  Poisson  breakdown  process  in  more  detail.  We  also  state 
conditions  under  which  the  numerical  procedures  developed  in 
Chapter  2 are  applicable  to  the  present  model.  Finally,  in 
section  3 . 5 we  suggest  two  approximation  methods. 
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3.2.  Derivation  of  the  Transform  of  the  Distribution 
of  Virtual  Service  Period 

Treated  as  a function  of  the  parameter  t,  the  virtual 
service  time  V(t)  is  a stochastic  process  that  represents 
the  total  time  required  to  achieve  t units  of  active  ser- 
vice. We  have  that 

V(t)  = t + Q(t)  (3-2-1) 

where  Q(t)  is  a process  called  excess  time,  introduced  by 
Muth  [1970,  1979a]. 

With  S the  service  time  of  a station,  the  virtual 
service  period  of  that  station  becomes 

V(S)  = S + Q(S)  (3-2-2) 

where  Q(S)  represents  the  accumulated  time  spent  on  repair 
during  that  service  period  (see  Figure  3.1). 

The  excess  time  Q(t)  is  a random  sum  of  breakdown 
periods.  We  define  N(t)  to  be  the  number  of  breakdowns  in 
the  interval  (0,t).  Its  generating  function  is 

CO 

G(t,z)  ^ I p(x,n)  (3-2-3) 

n=0 

where  p(x,n)  A P[N(x)=n]. 
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b 


Figure  3.1.  A particular  realization  of  Q(S)  and  V(S). 
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We  let  Tj  to  be  the  repair  time  following  the  j-th  break- 
down. We  assume  that  3.re  independently 

and  identically  distributed  random  variables,  each  of  them 
independent  of  N(t).  It  then  follows  that 

N(t) 

Q(t)  = E T (3-2-4) 

n=0 

where  ^ 0.  Using  conditioning  on  N(t)  the  distribution 
function  of  Q(t)  is  obtained  from  (3-2-4)  as 

00 

F (X)  = E F^^^Nx)  p(T,n)  (3-2-5) 

' n=0 

where  F^^^  (x)  is  the  n-fold  convolution  of  F^(x),  and 

(0)  * 

^ X i 0.  Let  F (s)  be  the  Laplace-Stielt jes 

transform  of  F(t),  then 


F 

Q(t) 


(s) 


E [e 


-sQ(t) j 


E [F*(s)]”  p(x,n) 
n=0 


(3-2-6) 


or,  using  (3-2-3) 


F*^^j  (s)  = G(t,  F*(s)  ) 


(3-2-7) 
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To  simplify  notation  we  now  let  V A V(S) . The  Laplace- 

* 

Stieltjes  transform  F^(s)  of  the  distribution  function  of  V 
is  obtained  in  two  steps.  The  first  step  leads  to  the 
transform  of  the  distribution  of  the  process  V(t).  From 
(3-2-1)  and  (3-2-6)  we  see  that 


* 

F 

V(t) 


(s) 


E[e"®(T+Q(T))] 


e"®^  G(t,F*(s) ) . 


In  the  second  step  we  condition  on  S = t to  obtain  the 
distribution  of  V.  Since 


F^(x)  = / (X)  dF  (t) 

T = 0 ' 


it  follows  that 


F*(s)  = e"®^  G(t,F* (s) )dFg(T)  (3-2-8) 

This  is  a completely  general  formula  relating  the  distri- 
bution of  the  virtual  service  period  to  the  distribution  of 
the  service  time,  the  distribution  of  the  repair  time  and 
the  distribution  of  the  breaJcdown  counting  process.  In 
order  to  reduce  the  integral  in  (3-2-8)  to  an  algebraic 
expression,  specific  assumptions  need  to  be  made  about  the 
distributions  of  these  quantities.  Moreover,  in  order  to 
compute  the  throughput  rate  using  the  numerical  procedures 
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developed  in  Chapter  2,  we  must  require  that  F^(t)  be  of  a 
special  phase  type  distribution.  The  question  therefore 
arises  as  to  what  restrictions  must  be  placed  on  the  dis- 
tributions of  S,  T and  N(x)  in  (3-2-8)  in  order  to  achieve 
this  objective.  In  the  process  of  specializing  (3-2-8)  we 
may  begin  by  first  selecting  a specific  distribution  for  S 
or  by  first  selecting  a specific  model  for  the  breakdown 
process  N(t).  The  latter  case  will  be  treated  in  section 
3.3. 

3.2.1.  The  Case  Where  Service  Time  is  of  Special  Phase  Type 
In  the  less  complex  model  of  Chapter  2 the  service 
distribution  was  restricted  to  that  of  special  phase  type, 
for  the  reason  just  mentioned.  It  is  therefore  appropriate 
to  apply  the  same  restriction  to  this  model.  We  then  have 


m 


1 - E 
i=l 


and  equation  (3-2-8)  becomes 


★ m 

F (s)  = E a. B.  / G(t 

i=l  0 


i=l 


E G (s  + a^,  (s)  ) . 


(3-2-9) 
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where  G(s,z)  is  the  Laplace  transform  of  G(t,z).  This 
equation  will  be  further  specialized  in  section  3.3  to  a 
model  of  N(x)  that  has  a closed  form  expression  for  G(s,z). 

3.2.2.  The  Simplest  Model 

As  an  example  of  how  (3-2-8)  may  be  reduced  to  an 
algebraic  expression  we  treat  the  simplest  model  which 
arises  when  S and  T are  exponentially  distributed  with  rates 
a and  3 and  when  N(x)  is  a Poisson  process  with  rate  A.  For 
this  case,  we  have 


Fs(t)  = 

F*(s)  = 
G(x,z)  = 


1-e 


-ax 


6 

s + 3 
g-Ax  (1-z) 


Substitution  into  (3-2-8)  yields 


Fv(s) 


a f 
0 


-sx 


e-A(l  - Ft(s))t  g-ax^^ 


a 

" - i-H' 

g (s  + 3) 

2 

s + (a  + 3 + A)s  + a3 


(3-2-10) 
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The  above  expression  was  previously  obtained  by  Altiok  and 

Stidham  [1983]  in  a direct  manner  by  using  the  no-memory 

property  of  exponential  random  variables. 

* 

The  poles  of  Fy(s)  can  be  obtained  from  equation 
(3-2-10)  as 


’1,2 


_ — (g+6+X)  f / (g+S+i)  ~ 4a 6 


(3-2-11) 


After  some  manipulation  we  have 


’1,2 


-(g+B+X)  ± /(g-B-A)  + 4aX 


(3-2-12) 


From  equations  (3-2-11)  and  (3-2-12)  we  can  see  that  s^  and 

* 

S2  are  always  real,  distinct  and  negative.  Therefore  F^(s) 
always  has  real,  distinct  and  negative  poles  for  all  combi- 
nation values  of  g,  B and  X.  This  means,  for  the  simplest 
model,  F^(t)  is  always  a special  phase  type  distribution. 

In  section  3.4  we  discuss  more  general  models  for  which 
F^(t)  is  still  a special  phase  type  distribution. 

3.3.  Model  for  N(t) 

In  this  section  we  select  specific  models  for  the 
breakdown  counting  process  and  substitute  them  into  the 
formulation  (3-2-8)  to  obtain  algebraic  expressions  of 
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X if 

F^(s) . We  also  state  the  condition  under  which  F^(s)  is  of 


phase  type. 

3.3.1.  Poisson  Process 

The  simplest  model  is  obtained  when  N(t)  is  a Poisson 
process.  For  this  model  we  have 


where  U is  the  time  between  breakdowns  and  \ is  the  station 
breakdown  rate.  Also 


P (t ,n) 


( At ) ^ -At 
n!  ® 


and 


G(t,2) 


e 


-At  (1-z) 


(3-3-1) 


Substitution  of  (3-3-1)  into  (3-2-8)  yields 


0 


= Fg(s  + A - A F^(s)) 


(3-3-2) 
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The  same  expression  was  obtained  by  Gaver  [1962]  as  the 
total  service  for  a customer  with  low  priority  in  an  M/G/1 
queueing  system  with  two  priorities.  The  repair  time  T in 
(3-3-2)  corresponds  to  the  service  time  of  a customer  with 
high  priority. 

3.3.2.  Renewal  Counting  Process 

If  the  sequence  of  time  between  breakdowns  , 

U2 , . . .is  independently  and  identically  distributed,  then 
this  sequence  forms  an  ordinary  renewal  process  and  N(t)  is 
the  corresponding  renewal  counting  process.  An  ordinary 
renewal  process  in  the  context  of  our  model  implies  that  the 
beginning  of  a virtual  service  period  coincides  with  a 
renewal  point,  which  means  that  a repair  has  just  been 
completed.  In  reality  the  beginning  of  a virtual  service 
period  will  fall  randomly  between  renewal  points.  Hence  an 
equilibrium  renewal  process  is  a more  appropriate  model. 

In  an  equilibrium  renewal  process  U2 , , . . . are 

independently  and  identically  distributed  as  U,  and  is 
the  forward  recurrence  time  of  the  process  with  distribution 
function 


F 


U 


1 


(x) 


_1 

E[U] 


T 

f (1  - Fy(x))dx 


(3-3-3) 


The  probability  law  of  N(t)  is  linked  to  the  proba- 
bility law  of  U and  by  the  well  known  relation 


64 


p(T,n)  = 


1 - F„^(x) 


for  n=0 


Fu^(x) 


* 

U 


(t)  * (1  - 


Fu  (T)) 


for  n=l ,2 , . . . 
(3-3-4) 


where  * denotes  convolution. 

Taking  the  Laplace  transform  of  (3-3-4)  gives 

(1  - Fy  (s)  ) for  n=0 

^F*^(s)  (F*(s))^"^  (1-F*(s))  for  n=l,2,... 

(3-3-5) 

Applying  the  definition  (3-2-3)  to  (3-3-5)  we  obtain  the 
transform  of  the  generating  function  of  N(x)  as 


p(s,n)  = 


zF*  (s)  (1-F*(s)) 


G ( s , z ) — — 


U, 


U 


(1-Fy  (s)  + 


1-z  Fy(s) 


or 


1-F*^ (s)  + z (F*^ (s)  - F*(s) ) 
s (1-z  F*  (s) ) 


G (s,z) 


(3-3-6) 
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where  from  (3-3-3) 


U. 


(s)  = 


l-Fy(S) 

“eTuI 


Thus,  after  some  manipulation  we  have 


T (1-z)  (1-F..(s)  ) 

G (s,z)  = y-j (3-3-7) 

sE[U] (1-z  F (s) ) 


If  is  identical  to  U then  we  have  an  ordinary  renewal 
process  and  (3-3-6)  reduces  to 


G (s,z) 


l-Fy(s) 


s (1-zFy (s) ) 


Substitution  of  (3-3-7)  into  (3-2-9)  yields 


* m (1-F*(s))  (l-F*(s+a  )) 

F (s)  = Z cx  3.  {— * ^ ^ ) (3-3-8) 

1=1  SE[U](1-F  (s)  F (s+a.)) 


From  this  equation  we  can  see  that  if  Fy(s)  and  F^(s)  are 

* 

rational  functions,  then  F^(s)  is  also  a rational  function. 
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3.3.3.  Phase-Renewal  Counting  Process 

This  special  case  of  a renewal  process  is  obtained  by 
letting  the  time  between  breakdowns  U be  a phase  type  random 
variable.  We  again  consider  the  case  of  an  equilibrium  pro- 
cess. For  simplicity,  assume  that  U is  a continuous  random 
variable.  In  order  to  make  subsequent  results  notationally 
compact,  the  distribution  function  of  U is  expressed  as 

Fy(T)  = 1 - q e'^'^u  (3-3-9) 

where  q is  a row  probability  vector,  u is  a column  vector 
whose  elements  are  1,  and  C is  a matrix  with  special  prop- 
erties as  discussed  in  appendix  A.  The  Laplace-Stielt j es 
transform  of  F^(t)  is 

F*(s)  = q (-C)  (sI-C)"^u  (3-3-10) 


From  (3-3-3)  and  (3-3-9)  we  obtain 


(s)  = 


q (sI-C)  ^u 
E[U] 


(3-3-11) 


By  comparing  (3-3-10)  with  (3-3-11)  it  is  apparent  that 

is  also  of  phase  type  and  that  F (t)  consists  of  a 

1 

combination  of  exactly  the  same  exponential  functions  that 
are  present  in  Fy(x).  Therefore  we  must  be  able  to  write 
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F*  (s)  = TT  (-C)  (sI-C)”^u 


(3-3-12) 


and 


U. 


(t)  = 1-TT 


Ct 

e u 


(3-3-13) 


where  ir  is  a row  probability  vector.  Comparing  (3-3-11)  and 
(3-3-12)  we  see  that 


TT  (-C) 


E[U] 


Multiplying  the  above  equation  by  the  matrix  uq  yields 


IT  (-C)  uq 


= 

E[U] 


E[U] 


and  therefore 


TT  (-C)  = TT  (-C)  uq. 


The  result  is  the  defining  equations  for  the  vector  tt  , 
namely 


IT  (C-Cuq) 


0 


IT  U 


1. 
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This  equation  for  (t)  was  derived  by  Neuts  [1981]  using 
a different  method. 

Substitution  of  Fy^(r)  and  Fy(x)  into  (3-3-4)  and  taking 
the  Laplace  transform  of  the  resulting  expression  yields 

( ir  (sI-C)  u for  n=0 

p(s,n)  = i 

^ TT  (Sl-C)"^(-C)u[q(sl-C)"^(-C)u]'""^q(sl-C)"^u 

for  n=l ,2. . . 


Writing  B = -Cuq,  this  equation  becomes 

p(s,n)  = TT  ( (sI-C)“^B)'^(sI-C)“^u  for  n 

Substitution  in  (3-2-3)  yields 

G (s,z)  = TT  (I-z  (sI-C)  ”^B)  (sI-C)  ~^u 
Since  B and  (sI-C)  commute,  we  have 

G (s,z)  = IT  (sI-C-zB)  ^u 

= ir  (sI-(C+zB))  ^u 


0,1,2. . . 
(3-3-14) 


(3-3-15) 


The  generating  function  of  N(t)  is  now  obtained  by  taking 
the  inverse  Laplace  transform  of  (3-3-15) 
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(3-3-16) 


or 


If  we  wanted  to  use  an  ordinary  renewal  process  we  would 
simply  replace  ir  in  the  above  equation  with  q.  Equation 
(3-3-16)  was  also  derived  by  Neuts  [1981]  using  a different 
method. 

As  a final  step  we  determine  the  transform  of  the 
distribution  of  virtual  service  period.  Substitution  of 
(3-3-16)  into  (3-2-8)  yields 


u dFg  (t) 


0 


dFg(x)  u 


0 


= TT  Fg  (sI-C-F^  (s)  B)  u 


(3-3-17) 


It  can  be  seen  from  (3-3-17)  that  F*(s)  will  be  a rational 


function  if  F^(s)  and  Fg(s)  are  rational  functions. 
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3.4.  The  Case  of  Poisson  Breakdown  Process 
It  is  reasonable  to  assume  that  the  station  time 
between  breakdowns  has  no-memory,  which  means  that  the 
breakdown  counting  process  is  a Poisson  process.  For  this 
case,  we  showed  that 


F*(s)  = F*(s  + A - A F*(s))  (3-3-2) 

ic 

As  discussed  in  section  3.2  we  wish  to  have  F^(s)  of  special 

phase  type  so  that  the  procedures  for  calculating  the 

throughput  rate  developed  in  the  previous  chapter  can  be 

applied.  In  this  section  we  determine  those  special  models 
* * * 

of  Fg(s)  and  F^(s)  that  lead  to  F^(s)  of  special  phase  type. 
From  the  structure  of  (3-3-2)  we  conclude  first  that 

* 

Fg(s)  needs  to  be  a rational  function  with  distinct  poles 
only.  Thus,  writing 


F*(s) 


m 

I 

i=l 


°-i^i 

(s+ai) 


we  obtain 


Fv(s)  = 


m 

Z 


a.,. 


s+a.+A-A  F„(s) 

1 ip  ' ' 


i=l 


(3-4-1) 
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Next  we  conclude  that  if  F^(s)  is  a rational  function  then 
* 

F^(s)  becomes  a rational  function.  However  the  poles  of 
* 

F^(s)  are  not  necessarily  all  real  even  though  all  a.  are 
★ 

real  and  F^(s)  has  only  real  poles.  It  is  too  difficult  to 

investigate  the  behavior  of  the  poles  of  F* *(s)  for  a 

general  model  of  F^{t).  Instead  we  will  investigate  the 
* 

poles  of  F^(s)  for  two  specific  models  of  F^{t),  one  is  a 
hyperexponential  distribution,  the  other  is  a two-stage 
general  Erlang  distribution,  that  is,  a hypoexponential 
distribution.  These  two  distributions  are  the  simplest  form 
of  special  phase  type  distributions  that  permit  the  value  of 
the  squared  coefficient  of  variation  to  be  different  from  1. 
The  coefficient  of  variation  of  a hyperexponential  dis- 
tribution is  greater  than  1 and  that  of  a two-stage  general 
Erlang  distribution  is  between  /O . 5 and  1. 


Case  1:  F^(t)  is  a hyperexponential  distribution 

For  this  case  we  have  the  following  result: 

Result  3.1;  If  T is  a hyperexponential  random  variable, 

•k 

then  all  poles  of  F^(s)  are  real,  distinct  and 
negative . 

Proof ; If  T is  a hyperexponential  random  variable,  then 

* 

F^(s)  has  the  form 


F,p(s)  = 


Y . 6 . 
1 1 


(s+6^) 


n 

Z 

i=l 


(3-4-2) 
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where  n is  finite; 


"i' 


6 . >0  and 
1 


n 

E 

i=l 


1. 


For  convenience,  let  the  6^  be  assigned  such  that 
<5]^  < 52  '^  * ’ * *^  Consider  first  the  j-th 

term  of  the  summation  in  equation  (3-4-1) , which  is 


h.  (s) 


a .8  . 

-J-J 


n Y • 5 . 

(s+a.  + A-A  E ' 

J (s+6.) 


9j  (s) 


The  denominator  of  this  expression  is  an  increasing 
function  of  real  s for  < s < —6.  (see  figure 

3.2)  and  therefore  has  exactly  one  real  zero  in  each 

such  interval.  It  follows  from  a similar  argument 
that  there  is  exactly  one  real  zero  in  the  interval 

(-"/  exactly  one  real  zero  in  the  interval 

(-6^,  0).  Since  gj  (s)  is  of  degree  (n+1)  there  can- 
not be  any  other  zeros.  Consider  next  the  zeros  of 
gj^(s),  where  k j . The  function  gj^(s)  differs  from 

g.  (s)  because  of  the  parameter  ^ , Thus  if  ^ 5^  2 _ 

3 k k j 

then  for  each  of  the  n+1  nonoverlapping  intervals  on 

the  s-axis,  the  zero  of  gj^(s)  is  different  from  the 

zero  of  g^ (s) . This  case  is  shown  in  figure  3.2. 

Hence  if  the  j = 2^  2,  ...,  m are  distinct,  then 

★ 

the  poles  of  F^(s)  are  distinct. 
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Figure  3.2. 


The  graphs  of  g.  (s)  and  g,  (s) 
J 

for  -<S . , - < s < -6  . . 

1+1  1 
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This  result  guarantees  that  the  numerical  procedures  for 
calculating  the  throughput  rate  developed  in  the  previous 
chapter  can  be  used  as  long  as  (t)  is  of  special  phase 
type  and  F^(t)  is  hyperexponential.  Note  that  if  we  wish  to 
prescribe  only  the  first  two  moments  of  the  repair  time, 
then  we  need  only  two  terms  in  the  sum  (3-4-2) . 

Case  2;  F^(t)  is  a two-stage  general  Erlang  distribution 

If  T is  two-stage  general  Erlang  random  variable,  then 


Consider  first  the  j-th  term  of  the  summation  in  equation 
(3-4-1) , which  is 


F^(s)  has  the  form 


(s+6^) (s+62) 


s+cx  . +A  — A 
3 


(s+6^) (s+62) 


The  poles  of  h^  (s)  are  the  roots  of  the  characteristic  equa- 
tion 


s (s+6^+62) 


(s+Oj) (s+6^) (s+62) 


1 + A 


0 


(3-4-3) 
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We  wish  to  determine  the  range  values  of  X,  a.,  6.  and  6 
for  which  these  poles  are  real.  A convenient  tool  for 
investigating  the  roots  of  an  equation  is  the  root  locus 
plot,  A brief  explanation  of  the  root  locus  method  follows. 
Let  equation  (3-4-3)  be  written  as 

1 + X G(s)  = 0 

where  X may  take  values  in  (0,  ~) . The  roots  of  this  equa- 
tion are  the  values  of  s that  satisfy 

s = g"^ (-1/X) 

The  root  locus  is  a one  dimensional  curve  that  maps  the 
negative  real  axis  -1/X  into  the  s-plane  via  the  mapping 
function  G and  parameterized  by  the  values  of  X.  Of 

interest  here  are  those  sections  of  the  root  locus  that  lie 
on  the  real  s-axis.  A reference  on  the  methodology  of  root 
locus  plotting  is  Dorf  [1967]  . 

In  applying  the  root  locus  method  to  equation  (3-4-3) 
it  turns  out  that  two  separate  cases  need  to  be  considered 
here.  In  both  cases,  the  root  locus  has  three  distinct 
branches  corresponding  to  the  three  solutions  of  (3-4-3) . 
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Case  1; 


“i  ^ ^ ^2 


The  root  locus  of  this  case  is  shown  in  figure  3. 3. a. 

★ 

It  IS  seen  that  the  poles  of  F^(s)  are  always  real,  distinct 
•and  negative. 


Case  2; 


a . 

_J_ 


6i  + 


The  root  locus  of  this  case  is  shown  in  figure  3.3.b. 
The  figure  shows  the  case  a.  >6.  If  6^  < a.  <6,  or 

3 ^ I j 2 

< 6^  the  structure  of  the  root  locus  plot  remains 
unchanged;  and  only  the  order  of  a.,  6 and  6 on  the  real 
axis  is  changed.  In  this  case  there  is  a range  of  values  of 
A in  which  two  of  the  poles  of  F^(s)  are  complex  conjugate. 
The  range  of  values  of  A that  still  give  real  poles  is 
obtained  as  follows.  From  equation  (3-4-3)  represent  A as  a 
function  of  real  s 


A = A (s) 


and  then  find  the  real  solutions  of  the  equation 


dA  (s) 
ds 


(3-4-4) 


In  our  case  the  above  equation  has  four  roots,  two  of  which 
are  real  and  the  other  two  are  complex  conjugate.  While  it 
is  easy  to  find  these  two  real  roots  numerically,  their 
closed  form  expressions  are  not  so  simple.  Let  the  real 
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X:  7\  = 0 


O : 7\=  oo 


Figure  3.3.  The  root  locus  of  equation  (3-4-3)  with  X as 
the  parameter. 
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solution  of  (3-4-4)  be  and  S2  and  X(s^)  < X(s2),  then  the 
range  values  of  X that  gives  real  roots  are 

0 S X ^ X(s^)  and  X ^ ^^^2^ 

The  following  result  is  a consequence  of  the  first  case. 
Result  3.2;  If  T is  a two-stage  general  Erlang  random 

variable  and  min  {a.}  ^ 6 + & , then  the  poles 

* J 

of  F^(s)  are  all  real,  distinct  and  negative. 
From  the  above  discussion  we  see  that  the  poles  of  F*(s)  are 
not  always  real  when  F^(t)  is  a two-stage  Erlang  distribu- 
tion. Obviously  the  possiblity  of  obtaining  complex  poles 
will  be  greater  for  more  complex  models  of  F^(t),  such  as 
three-  or  four-stage  Erlang  distributions. 

3.5.  Approximation  Methods 
From  equation  (3-3-17)  we  can  see  that  the  degree  of 
complexity  of  the  virtual  service  period  distribution  is  the 
compounded  effect  of  the  complexities  of  the  distribution 
of  the  service  time,  the  repair  time  and  the  time  between 
breakdowns.  In  practice  these  random  variables  do  not  have 
a simple  distribution  function,  such  as  exponential,  since 
their  coefficients  of  variations  are  usually  very  different 
from  1.  The  transform  of  these  distributions  therefore  may 
have  a large  number  of  poles.  If  this  is  the  case  then  the 
virtual  service  period  distribution  becomes  very  complex, 
meaning  that  the  poles  of  its  transform  are  very  likely  to 
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be  nonreal  or  nondistinct.  Even  when  it  is  reasonable  to 
assume  that  the  time  between  breakdowns  is  exponential,  the 
service  time  and  the  repair  time  may  still  have  complex 
distribution  functions,  and  it  was  shown  that  for  this 
simple  model  of  N(t)  the  poles  of  the  transform  of  the 
virtual  service  period  distribution  are  not  always  real  and 
distinct.  Even  if  all  poles  are  real  and  distinct  there  may 
be  so  many  of  them  as  to  make  the  numerical  procedure 
inefficient. 

Another  difficulty  is  that  in  complex  cases  the  primary 
information  we  have  about  the  random  variables  service  time, 
repair  time  and  time  between  breakdowns  usually  is  not  on 
their  distribution  functions,  but  only  limited  on  their 
first  two  moments. 

To  approach  such  complex  cases,  it  is  probably  better 
to  use  efficient,  but  fairly  accurate,  approximation 
methods,  rather  than  to  insist  on  exact  but  cumbersome 
computations.  Two  such  approximation  methods  are  presented 
next.  One  is  the  "single  breakdown  approximation"  and  the 
other  is  the  "two  moment  approximation." 

3.5.1.  Single  Breakdown  Approximation 

In  order  for  a system  to  be  usable  in  practice,  the 
probability  of  one  or  more  breakdowns  during  any  single 
service  period  should  be  small.  As  a consequence,  the  prob- 
^t>ility  of  more  than  one  breakdown  during  any  service  period 
will  be  negligibly  small.  Therefore,  it  is  reasonable  to 
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assume  that  during  a single  service  period,  a station  will 
break  down  at  most  once.  Such  an  approximation  was  sug- 
gested by  Muth  [1979a] . 

Let  be  the  time  until  the  first  breakdown  with 

distribution  function  F (t) . It  follows  that  the  breakdown 

^1 

counting  process  under  the  single  breakdown  approximation 
has  the  form 

(1  - F (t)  for  n=0 

Fy  (t)  for  n=l 

0 otherwise  (3-5-1) 

This  approximation  is  good  only  for  the  values  of  x such 
that  P[N(x)  = 0]  is  close  to  one.  As  a rule  of  thumb  we 
can  pick  a value  x for  which  P[N(x  ) = 0]  = 0.95  and  say 
that  the  approximation  is  justified  when  E[S]  + 2/Var (s)  < x*. 

The  generating  function  of  N(x)  is  obtained  from 
(3-5-1)  as 

G(x,z)  = 1 - F (x)  + F (x)z  (3-5-2) 


If  has  the  phase  type  distribution  given  under  (3-3-13) , 
then 


= ir 


Cx 

e u 


+ Z (l-TT 


Cx  . 
e u)  . 


G ( X , z) 


(3-5-3) 
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Substitution  of  the  above  expression  into  (3-2-8)  yields 


F (s)  = / e ir  e*"^u  dF„(T) 

0 


+ 


F (s)  / e (1-ir  e*"^u)  dF  (t) 

0 ^ 


After  some  simplification  we  obtain 

F*(s)  = F*(s)  F*(s)  + (1-F*(s))  IT  F*(sl-C)u  (3-5-4) 

Close  inspection  of  the  above  equation  shows  that  the  number 

of  poles  is  not  fixed  but  is  less  than  or  equal  to 

m(n  + 1)  + k where  m and  k are  the  number  of  poles  of  F„(s) 

D 

★ 

and  F^(s)  and  n is  the  size  of  matrix  C.  By  comparison,  the 

number  of  poles  resulting  from  the  exact  formula  (3-3-17) 

is  mn(k  + 1) . In  most  cases  mn(k  +1)  is  greater  than 

★ 

m(n  +1)  + k.  Let  Fg(s)  be  a special  phase  type  distri- 
bution with  poles  -a^,  i=l,  . . . , m and  F*(s)  be  a special 

phase  type  distribution  with  poles  -y . , j=l,  . . k.  We 

* 

then  have  that  the  poles  of  F^(s)  are  real  and  distinct  if 
and  only  if  and  y j are  distinct  for  all  i and  j , and  the 
n roots  of  the  equation  | (s+a^)I  - c|  = 0 are  all  real  and 
distinct  for  all  i and  not  equal  to  or  y^.  We  have  the 
following  result. 
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Result  3.3;  If  T and  S are  special  phase  type  random 

variables  with  poles  -a . and  -y  • where  a . ^ y . 

1 D I D 

for  all  i and  j , and  the  breakdown  process 

is  a Poisson  process  with  rate  X such  that 
^ Yj  for  all  i and  j,  then  the  poles 

★ 

of  F^(s)  are  all  real  and  distinct. 

Proof;  If  the  breakdown  counting  process  is  a Poisson 
process  with  rate  X,  equation  (3-5-4)  becomes 

F*(s)  = F*(s)  F*(s)  + (1-F*(s))  F*(s+X) 
and  the  proof  is  at  hand. 

Thus  to  attain  distributions  with  real  and  distinct 
poles  is  easier  with  the  approximation  (3-5-4)  than  with  the 
exact  formula  (3-3-17) . 

In  some  cases,  however,  the  exact  formula  (3-3-17)  is 
to  be  preferred  over  the  approximation  (3-5-4) . As  an 
example,  suppose  that  the  breakdown  counting  process  is  a 
Poisson  process,  the  time  to  repair  is  exponential  and  the 
service  time  is  of  special  phase  type.  If  the  exact  formula 
is  used  then  the  number  of  poles  will  be  2m  and  according  to 
result  3.1  all  of  them  are  real,  distinct  and  negative.  If 
the  approximation  formula  is  used  then  the  number  of  poles 
will  be  2m  + 1 and  all  of  them  are  real  but  not  necessarily 
distinct.  For  the  poles  to  be  distinct,  using  result  3.3, 
the  repair  rate  minus  the  breakdown  rate  must  be  not  equal 
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★ 

to  any  poles  of  Fg(s).  This  approximation  method  is  useful 
only  if  the  exact  method  fails,  because  it  does  not  yield 
real  poles,  or  if  mn(k+l)  > m(n+l)  + k and  the  approximation 
gives  real  and  distinct  poles.  If  the  exact  and  the  single 
breakdown  approximation  methods  fail,  then  the  two  moment 
approximation,  which  will  be  given  in  the  next  section, 
should  be  used. 

3.4.2.  Two  Moment  Approximation 

Muth  [1977]  computes  and  compares  the  throughput  rate 
of  three-station  balanced  production  lines  when  all  stations 
have  Erlang  service  times  and  when  all  stations  have  uniform 
service  times.  The  comparison  is  made  with  distributions 
whose  first  two  moments  are  matched.  The  results  show  that 
the  throughput  rates  over  a range  of  values  of  the  coeffi- 
cient of  variations  are,  for  all  practical  purposes,  iden- 
tical. Since  these  two  distributions  differ  considerably  in 
their  higher  moments,  the  results  indicate  that  the  third 
and  higher  moments  of  the  service  time  have  negligible 
effect  on  the  throughput  rate.  Applying  this  result  to  our 
case  we  conclude  that  the  throughput  rate  of  a production 
line  subject  to  breakdown  may  be  approximated  by  a method 
that  requires  as  input  only  the  first  two  moments  of  the 
virtual  service  periods. 

Such  an  approximation  is  useful  if  the  virtual  service 
period  V has  a distribution  function  that  is  either  very 
complex,  meaning  the  transform  has  a large  number  of  poles. 
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or  if  the  transform  has  complex  or  nondistinct  poles.  The 
key  idea  of  the  approximation  is  to  compute  the  first  two 
moments  of  V directly  from  the  parameters  of  S,  T and  N(x), 
and  then  to  compute  the  throughput  rate,  using  the  method- 
ology of  chapter  2,  with  the  simplest  and  the  most  conven- 
ient distribution  of  service  times  that  matches  the  first 
two  moments  of  V. 

First  we  need  to  establish  general  expressions  that 
relate  the  first  two  moments  of  V to  the  parameters  of  S,  T 
and  N(t).  One  way  to  obtain  these  expressions  is  by  using 
the  formula 


where  F^(s)  is  given  by  equation  (3-2-8).  Here  we  use  a 
more  direct  approach  as  in  Muth  [1979a] . We  have 


E[V"^]  = (-1)^ 


* 


V = S + Q(S) 


(3-2-2) 


or 


N(S) 

V = S + I T. 

• 1 1 


(3-5-5) 
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Since  the  are  independently  and  identically  distributed 
as  T,  taking  the  expected  value  of  the  two  sides  of  equation 
(3-5-5)  yields 

00 

E[V]  = E[S]  + Z i E[T]  P[N(S)  = i] 
i=l 


or 


E[V]  = E[S]  + E[T]  E[N(S)  ] 

By  conditioning  E[N(S)]  on  S we  obtain 

00 

E[V]  = E[S]  + E[T]  f E[N(t)]  dF„(T)  (3-5-6) 

0 ^ 

2 

The  second  moment  E[V  ] is  obtained  by  squaring  equation 
(3-2-2) , completing  the  square  on  the  right  hand  side,  and 
taking  the  expected  value  of  all  terms.  The  result  is 

E[V^]  = E[S^]  + 2E[SQ(S)]  + E[Q(S)^] 


and,  conditioning  on  S as  before,  we  obtain 
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E[V^]  = E[S^]  + 2E[T]  / T 

0 


E[N(t)]  dFg(T)  + 


00 

Var[T]  / E[N(t)]  dF„(T) 

0 ^ 


+ E[T]^  / E[N(t)^]  dF^(x) 

0 ^ 


(3-5-7) 


If  N(t)  is  a Poisson  process,  then  we  have  the  following 
formulas : 


= Fbr 

e[n(t)  — 

E[U]  E[U] 

Substitution  of  the  above  equations  into  (3-5-6)  and  (3-5-7) 
yields,  after  some  simplifications, 

E[V]  = E[S]  (1  + E[T]/E[U])  (3-5-8) 

E[V^]  = E[S^]  (1  + E[T]/E[U])^ 


+ (1  + n^) 


E[T]  E[S]^ 
E[U] 


(3-5-9) 


where  n denotes  coefficient  of  variation. 

From  these  two  equations  we  can  see  that  we  need  to 
have  information  only  about  the  first  two  moments  of  S,  T, 
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and  U in  order  to  compute  the  first  two  moments  of  V.  Here 
we  are  interested  in  the  dependence  of  n.,  on  This 

relation  will  let  us  select  an  appropriate  distribution  that 
matches  the  first  two  moments  of  V.  After  some  simplifica- 
tion we  obtain  from  (3-5-8)  and  (3-5-9)  the  following 
relation : 

ny  = Hg  + (1  + n^)  g(x,y)  (3-5-10) 

where  2 

g(x,y)  = - ^ - 

(x+y)^ 

X = E [T]  /E  [S] 
y = E[U]  /E[S] 

2 2 
It  IS  seen  that  Hy  is  always  greater  than  or  equal  to  rig. 

Figure  3.4  shows  the  contours  of  g(x,y)  as  x and  y are 

varied  from  0 to  5.0.  Using  these  contours  riy  can  be  found 

easily. 

Now  we  select  the  simplest  model  for  the  moment- 
matching distribution  to  be  used  in  the  calculation  of  the 
throughput  rate.  Depending  on  the  value  of  the  squared 
coefficient  of  variation  of  the  virtual  service  period,  two 
models  of  special  phase  type  distribution  functions  can  be 
selected.  One  is  the  hyperexponential  distribution,  for  the 
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Figure  3.4.  Contours  of  constant  g(x,y). 
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case  > 1,  and  the  other  is  the  hypoexponential,  for  the 
2 

case  riy  2 1. 


Case  1;  The  virtual  service  period  is  hyperexponential 
The  simplest  model  is  a mix  of  only  two  exponential 
functions,  namely 


F^(t)  = 1 - 3e 


-ait 


- (1-6)  e 


This  model  is  sufficient  to  cover  the  range  of  all  values  of 
2 

Hy  greater  than  1.  We  must  determine  the  parameters  a^,  a2 
and  6 that  satisfy  the  equations 

^ ^ (3-5-11) 

1 2 

^ = E[V^]/2  (3-5-12) 

“l  “2 

^1'  “l'  “2  ^ ° 

These  two  equations  in  three  unknowns  do  not  have  a unique 
solution.  In  practice  one  may  choose  one  parameter  and  then 
determine  the  other  two  from  the  above  two  equations.  One 
way  to  do  this  is  by  imposing  the  following  restriction: 


90 


i_  = Izi 

“l  “2 


(3-5-13) 


This  restriction  means  that  each  exponential  term 
contributes  exactly  one-half  of  the  mean  value.  From 
equations  (3-5-11)  through  (3-5-13),  we  obtain 


_ 26 
1 E[V] 

and 


a 


2 


a 


1 


Case  2;  The  virtual  service  period  is  a hypoexponential 

The  simplest  model  for  this  case  is  an  m-stage  general 
Erlang  distributions  with  the  smallest  number  of  stages 
required  to  attain  the  desired  coefficient  of  variation  riy. 
We  thus  have 


Fy(t) 


m -a  . t 

1 - E 6.  e ^ 
i=l  ^ 
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where 


m 

n 


m 


n 


j=i 

j^i 


and  m is  the  smallest  integer  greater  than  — i— 


'v 


The  parameters  must  satisfy  the  following  equations: 


m 

E 

i=l 


E[V] 


(3-5-14) 


m 

E 

i=l 


1 

a . 2 
1 


Var  [V] 


(3-5-15) 


> 0 for  i = 1 , 2 , . . . , m 


We  have  two  equations  with  m unknowns.  These  equations  do 
not  have  a unique  solution.  One  way  to  cooose  a suitable 
set  of  is  to  impose  the  following  restrictions: 


_1 1^  ^ aE  [V] 

“i+1  “i  ^ 


i = 1,  2, 


m-1  (3-5-16) 


where  a is  some  postive  constant  that  needs  to  be  determined 
so  that  equations  (3-5-14)  and  equation  (3-5-15)  are 


satisfied.  From  equations  (3-5-14)  through  (3-5-16)  we 
obtain 
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a = 


3 (mriy  - 


m 


1) 

1 


and 


I_  = EtV] 

a . m 
1 


(1  + a(i 


m + 1 


) ) 


CHAPTER  4 


THE  BOWL  PHENOMENON 
4.1.  Introduction 

The  problem  of  unbalancing  a production  line  subject 
to  certain  constraints,  which  we  refer  to  as  constrained 
unbalancing,  has  attracted  attention  in  almost  25  papers 
since  it  was  introduced  by  Hillier  and  Boling  [1966]. 

Hillier  and  Boling  analyzed  a three-station  production 
line  with  exponential  service  times,  in  which  two  stations 
are  identical,  subject  to  the  constraint  that  the  sum  of 
the  mean  service  times  of  the  three  stations  is  constant. 
They  find  that  the  maximum  throughput  rate  is  obtained  by 
decreasing  the  mean  service  time  of  the  middle  station  and 
by  increasing  the  mean  service  time  of  the  two  identical 
outside  stations.  Hillier  and  Boling  coined  the  term  "bowl 
phenomenon"  for  this  optimal  arrangement  of  mean  service 
times  along  the  line. 

Presumably  the  discovery  by  Hillier  and  Boling  has  the 
following  practical  implication.  If  the  total  amount  of 
work  to  be  performed  on  an  item  is  considered  a fixed 
quantity  that  can  be  allocated  to  the  individual  stations  at 
the  discretion  of  the  line  operator,  then  the  operator 
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should  deliberately  unbalance  the  line  by  assigning  to  the 
interior  station  a workload  that  is  less  than  one-third  of 
the  total  load.  This  has  the  effect  of  decreasing  the  mean 
service  time  at  the  interior  station.  The  unbalancing,  if 
done  in  the  right  amount,  will  increase  the  throughput  rate 
of  the  production  line. 

It  is  well  known  that  in  the  case  where  all  stations 
have  deterministic  service  times  the  station  with  the 
greatest  service  time  is  the  bottleneck  that  determines  the 
line  throughput  rate;  mean  value  unbalancing  in  this  case 
will  therefore  decrease  the  line  throughput  rate.  It  is 
also  known  that  the  line  throughput  rate  decreases  with 
increasing  service  time  variabilities,  i.e.  with  increasing 
second  moments  of  the  service  time  distributions.  Further- 
more, the  effect  of  service  time  variability  on  line  effi- 
ciency is  very  different  in  nature  from  the  effect  of 
unbalance  on  line  efficiency.  It  must  therefore  be  surmised 
that  the  bowl  phenomenon,  as  observed  by  Hillier  and  Boling, 
is  brought  about  firstly  by  the  strong  dispersion  of  the 
exponential  distribution  and  secondly  by  the  effect  of 
unbalancing  variances  which,  for  the  exponential  distri- 
bution, are  necessarily  coupled  to  the  means.  This  is 
contrary  to  the  conclusion  of  Hillier  and  Boling  which 
attributes  the  bowl  phenomenon  only  to  the  unbalance  of  the 
mean  service  times. 

In  order  to  shed  light  on  the  mechanisms  by  which  mean 
value  unbalance  and  variance  unbalance  affect  the  throughput 


95 


rate,  it  is  more  meaningful  to  have  a method  of  analysis  in 
which  the  mean  value  and  the  variance  of  each  station  can  be 
specified  separately.  One  particular  question  of  interest 
is  to  determine  the  effect  of  unbalancing  the  mean  service 
times,  with  the  variance  of  each  station  remaining  fixed, 
subject  to  the  constraint  that  the  sum  of  the  mean  service 
times  is  constant.  This  analysis  of  the  bowl  phenomenon  is 
now  possible  through  the  method  given  in  Chapter  2. 

The  chapter  is  organized  as  follows.  In  the  next 
section  we  review  some  of  the  literature  on  the  bowl  phe- 
nomenon that  is  relevant  to  our  work.  In  section  4.3  we 
discuss  several  models  for  constrained  unbalancing. 

Finally,  in  section  4.4,  we  present  numerical  studies  and 
the  conclusions.  Appendix  B contains  the  details  of  the 
methodology  of  Chapter  2 as  applied  specifically  to  the 
model  of  section  4.4. 


4.2.  Past  Work 

Since  the  paper  by  Hillier  and  Boling  in  1966,  there 
have  been  many  studies  devoted  to  the  investigation  of  the 
bowl  phenomenon.  However,  in  the  interest  of  brevity  we 
will  discuss  only  those  studies  that  relate  to  our  work. 

The  model  of  constrained  unbalancing  of  Hillier  and 
Boling  [1966]  is  based  on  a constant  sum  of  the  mean  service 
times,  with  the  coefficient  of  variation  of  service  times  at 
each  station  equal  to  1 . 
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Using  the  integral  formulation  for  the  throughput  rate 
of  three-station  production  lines  derived  by  Muth  [1973], 

Rao  [1976]  extended  the  analysis  of  the  bowl  phenomenon  to 
the  case  where  not  every  station  has  the  same  coefficient  of 
variation.  In  his  analysis  he  treats  three  different  models 
of  three-station  production  lines.  In  the  first  model,  two 
stations  have  exponential  service  times  and  one  has  fixed 
service  time.  In  the  second  model,  two  stations  have  fixed 
service  times  and  one  has  exponential  service  time.  In  the 
third  model,  the  two  outside  stations  have  exponential 
service  times  and  the  middle  station  has  uniform  service 
time.  He  constrains  the  sum  of  the  mean  service  times  to  be 
constant  and  the  coefficient  of  variation  at  each  station  to 
be  constant.  He  finds  from  his  numerical  results  that  the 
throughput  rate  can  be  increased  by  assigning  higher  mean 
service  times  to  the  stations  with  the  lower  coefficient  of 
variation  of  service  times.  He  calls  this  optimal  workload 
allocation  "variability  imbalance."  For  the  case  where  the 
two  outside  stations  have  exponential  service  times  and  the 
middle  stations  has  fixed  service  times,  this  implies  that 
the  line  throughput  rate  can  be  increased  by  assigning 
higher  mean  service  times  to  the  middle  station  than  to  the 
two  outside  stations.  The  optimal  arrangement  of  mean 
service  times  along  the  line  is  then  in  the  form  of  an 
inverted  bowl. 

In  some  cases  the  effect  of  variability  imbalance  and 
of  mean  value  unbalance  may  cancel  one  another  as  it  happens 


97 


when  the  two  outside  stations  have  exponential  service  times 
and  the  middle  station  has  uniform  service  time  with  a 
coefficient  of  variation  equal  to  0.5. 

Hillier  and  Boling  [1979]  describe  the  bowl  phenomenon 
for  symmetrical  production  lines  with  up  to  six  stations  and 
coefficient  of  variation  of  service  times  equal  to  l/Zin 
where  m=l,2,  . . . ,7.  They  assume  that  the  service 
times  of  all  stations  are  Erlang  with  the  same  coefficient 
of  variation.  This  choice  of  service  time  distribution  lead 
them  to  the  same  model  of  constrained  unbalancing  as  in 
their  earlier  paper  [1966] . In  all  cases  considered,  the 
line  throughput  rate  can  be  increased  by  assigning  mean 
service  times  that  decrease  toward  the  middle  of  the  line. 

El-Rayah  [1979a]  considers  production  lines  with  3,  4 
and  12  stations.  He  treats  two  different  cases  for  each 
model;  either  all  stations  have  normal  service  times,  or  all 
stations  have  lognormal  service  times.  He  constrains  the 
sum  of  the  mean  service  times  to  constant  and  the 
coefficient  of  variation  at  each  station  to  be  constant. 
Using  a simulation  method,  he  calculates  the  throughput  rate 
for  four  different  arrangements  of  mean  service  times:  a 

balanced  line;  a line  with  an  increasing  order  of  mean 
service  times;  a line  with  alternate,  low  and  high,  mean 
service  times;  and  a line  whose  middle  stations  having  a 
lower  mean  service  time  than  outside  stations.  Using  a 
statistical  test,  he  compares  the  simulation  results  of 
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these  four  arrangements  and  draws  the  conclusion  that  the 
greatest  throughput  rate  is  obtained  when  middle  stations 
have  a lower  mean  service  time  than  outside  stations. 

In  a second  paper,  El-Rayah  [1979b]  consider  the  same 
models  of  production  lines  as  in  his  first  paper  but  using 
different  constraints.  He  constrains  the  sum  of  the 
coefficient  of  variation  of  service  times  to  be  constant  and 
the  mean  service  time  at  each  station  to  be  constant  and 
equal  to  1.  Using  a simulation  method,  he  calculates  the 
throughput  rate  of  production  lines  for  four  different 
arrangements  of  coefficient  of  variations:  a balanced  line; 

a line  with  an  increasing  order  of  coefficient  of 
variations;  a line  with  alternate,  low  and  high,  coefficient 
of  variations;  and  a line  whose  middle  stations  having  a 
lower  coefficient  of  variation  than  outside  stations.  Using 
a statistical  test,  he  compares  the  simulation  result  of 
these  four  arrangements  and  draws  the  conclusion  that  the 
greatest  throughput  rate  is  obtained  when  middle  stations 
have  a lower  coefficient  of  variation  than  outside  stations. 

For  the  results  of  both  papers,  it  must  be  noted  that 
the  accuracy  attainable  by  simulation  methods  is  outside  the 
range  of  change  in  the  throughput  rate  that  is  known  from 
analytical  methods.  For  example,  in  Hillier  and  Boling 
investigation  the  maximum  percent  of  improvement  obtained  by 
unbalancing  a four-station  production  line  is  less  than  1%. 
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4.3.  Formulations  of  Constrained  Unbalancing 
From  the  result  of  Muth  [1977]  as  discussed  in  sec- 
tion 3.5.2,  we  know  that  the  system  throughput  rate  is 
affected  mainly  by  the  first  two  moments  of  the  service 
times  and  that  the  effect  of  higher  moments  is  negligible. 
It  is  therefore  appropriate  to  base  the  analysis  on  any 
suitable  two-parameter  distribution  of  service  times,  to 
select  the  two  parameters  in  accordance  with  the  desired 
mean  and  variance,  and  then  to  treat  the  resulting  through- 
put rate  as  if  it  were  a function  of  the  first  two  moments 
only;  that  is 


r = r(y,a) 

where  y is  a vector  whose  i-th  element  is  y^,  the  mean 
service  time  of  station  i and  a is  a vector  whose  i-th 
element  is  a . , the  standard  deviation  of  the  service  time  of 
station  i. 

The  problem  that  will  be  addressed  is  to  establish  the 
functional  dependence  of  r on  y and  a.  In  particular,  we 
wish  to  find  the  maximum  throughput  rate  for  the  case  where 
y and  a are  varied  independently  subject  to  certain  con- 
straints on  these  two  parameters. 

The  constraints  are  not  true  physical  constraints,  but 
are  rather  the  assumptions  of  a specific  model.  Therefore, 
several  formulations  are  possible.  For  example,  as  applied 
to  the  first  moment,  we  may  select  the  constraint  as 
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K 

Z 

i=l 


C 


1 


or 


K 

Z 

i=l 


1/m. 


C 


1 


where  K is  the  number  of  stations  on  the  line  and 
positive  constant.  Similarly,  applied  to  the  second 
we  may  select  the  constraint  as 


or 


or 


or 


C 


2 


C 


2 


K 

Z 1/a.  = C 
i=l  ^ 


E 1/a.  = C 

i=l  ^ 


where  C2  is  some  positive  constant. 


is  some 
moment , 
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We  can  formulate  several  models  of  constrained 
unbalancing  by  using  different  combinations  of  the  con- 
straints given  above.  In  the  numerical  study  of  section  4.4 
the  following  two  constraint  sets  are  used. 


K 

I 

i=l 


C 


1 


K 

I 

i=l 


a . 
1 


C 


2 


Pi, 


a . SO 
1 


(4-3-1) 


and 


K 

Z 1/y.  = 
i=l  ^ ^ 


K 

Z 1/a.  = C, 
i=l  ^ ^ 


> 0 (4-3-2) 

In  comparison,  the  model  treated  by  Hillier  and 
Boling  [1966,  1979],  Rao  [1976]  or  El-Rayah  [1979a]  has  the 


constraint  set 
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K 

E 

i=l 


^ 0,  ^ 0.  (4-3-3) 

In  Rao's  formulation  the  €2^  are  not  the  same  for  all  i;  in 
Hillier  and  Boling's  and  El-Rayah's  formulations  the  €2^  are 
the  same  for  all  i.  As  we  can  see,  in  formulation  (4-3-3)  p 
and  o are  not  varied  independently. 

The  model  of  Hillier  and  Boling  and  El-Rayah  is  for 

K 

C„ . = C„ . When  combined  with  the  constraint  E p . = , 

K 

the  constraint  = ^2^^  implies  that  E = '“2*“1  thus 

(4-3-3)  becomes  a special  case  of  the  constraint  set 
(4-3-1).  In  other  words,  the  model  of  Hillier  and  Boling 
and  El-Rayah  belong  to  the  class  of  models  defined  by 
(4-3-1)  with  the  added  restriction  that  p and  a are  coupled, 
instead  of  being  selectable  independently.  Rao  does  not 
K 

Jceep  E a ^ equal  to  a constant  when  he  unbalances  the  mean 

values.  As  a consequence,  his  results  are  in  contradiction 
to  those  of  Hillier  and  Boling  as  will  be  shown  in  section 


4.4. 
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To  investigate  the  effect  of  unbalancing  the  mean 
service  times  alone,  one  should  use  the  constraint  set 


K 

E 

i=l 


C 


1 


cr 


i 


C 


2i 


y . > 0 

1 


(4-3-4) 


Similarly,  to  investigate  the  effect  of  unbalancing  the 
service  time  variabilities  alone,  one  should  use  the 
constraint  set 


K 

Z 

i=l 


a . 


1 


C 


2 


^ 0,  i 0 (4-3-5) 

The  constraint  set  used  by  El-Rayah  [1979b]  is  the  same  as 
the  constraint  set  (4-3-5) . Note  that  both  (4-3-4)  and 
(4-3-5)  are  just  special  cases  of  (4-3-1) . 

It  is  reasonable  to  assume  that  the  throughput  rate  in 
the  above  models  has  only  one  maximum.  Hence,  using  the 
reversibility  property  of  production  lines  (Muth  [1979b]) 
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the  line  with  maximum  throughput  rate  will  be  a symmetric 
line  (see  Hillier  and  Boling  [1979]  for  a discussion).  A 
symmetric  line  is  a line  with  y . = y.,  . , , and  o . = . 

Since  we  are  interested  in  finding  a resource 
allocation  that  yields  the  maximum  throughput  rate,  in  the 
investigation  we  need  to  consider  only  the  case  where  the 
line  is  symmetrical.  By  doing  this,  we  cut  the  number  of 
decision  variables  by  more  than  a half  and  therefore  we 
reduce  the  computational  effort  considerably. 


4.4.  Numerical  Study 

In  this  numerical  study  we  treat  a three-station 
symmetrical  production  line  whose  stations  have  no  buffers 
and  are  not  subject  to  breakdown.  This  is  the  simplest  case 
for  which  the  bowl  phenomenon  occurs.  Station  1 and 
station  3 are  identical;  that  is  y^  = y^  and  = a^.  After 
imposing  constraints  on  the  sum  of  the  first  moment  and  the 
sum  of  the  second  moment  in  accordance  with  (4-3-1)  or 
(4-3-2)  we  are  left  with  two  decision  variables,  namely  the 
ratios  Uj^/U2 

The  throughput  rate  r is  obtained  using  the  methodology 
of  Chapter  2 applied  to  this  specific  three-station  case. 

An  expression  for  r is  given  as  equation  (B-14)  in 
Appendix  B.  Using  this  expression,  the  throughput  rate  r 
can  be  computed  for  arbitrary  distributions  of  the  service 
time  of  the  middle  station  and  for  special  phase  type 
distributions  of  the  service  times  of  the  two  end  stations. 


105 


A computer  program  that  uses  this  method  was  written. 
This  program  computes  the  throughput  rate  r over  a specified 
range  of  values  of  the  two  parameters  x = 

y = o^lo^-  Several  cases  are  investigated.  These  cases  are 
shown  in  the  tables  below. 


Table  4.1.  Case  l.a. 


Station  1 

station  2 

Station  3 

service  time 
distribution 

Erlang 

deterministic 

Erlang 

^i 

(balanced 

condition) 

1 

1 

1 

a . 
1 

(balanced 

condition) 

00 

• 

o 

0 

0.8 

Constraints:  The  sum  of  the  mean  values  is  constant. 

The  variance  of  each  station  is  constant. 
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Table  4.2.  Case  l.b.  (Rao's  Formulation) 


Station  1 

Station  2 Station  3 

service  time 
distribution 

exponential 

deterministic  exponential 

^i 

(balanced 

condition) 

1 

1 1 

a . 

1 

(balanced 

condition) 

1 

0 1 

Constraints : 

The  sum  of  the  : 
The  coefficient 
is  constant. 

mean  values  is  constant, 
of  variation  of  each  station 

Table  4.3.  i 

Case  2 . 

Station  1 

Station  2 

station  3 

service  time 
distribution 

Erlang 

Erlang 

Erlang 

‘^i 

(balanced 

condition) 

1 

1 

1 

o . 
1 

(balanced 

condition) 

0.8 

o 

• 

00 

0.8 

Constraints:  The  sum  of  the  mean  values  is  constant. 

The  sum  of  the  standard  deviations  is 
constant. 
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Table  4.4.  Case  3. 


Station  1 

Station  2 

Station  3 

service  time 
distribution 

Erlang 

Erlang 

Erlang 

"i 

(balanced 

condition) 

1 

1 

1 

o . 
1 

(balanced 

condition) 

0.8 

0.8 

0.8 

Constraints:  The  sum  of  the  reciprocal  mean  values  is 

constant. 

The  sum  of  the  reciprocal  standard  deviations 
is  constant. 


The  results  of  case  l.a  and  case  l.b  are  shown  in 
figure  4.1  where  the  relative  throughput  rate  is  plotted 
as  a function  of  the  ratio  The  relative  throughput 

rate,  r^^,  is  defined  as  r/rj^  where  rj^  is  the  throughput 
rate  of  the  balanced  line.  From  this  figure  we  can  see 
that  if  the  variance  of  each  station  is  kept  constant  while 
unbalancing  the  mean  service  times  (case  l.a)  then  the  bowl 
phenomenon  occurs,  i.e.  the  maximum  throughput  rate  is 
obtained  by  assigning  a lower  mean  service  time  to  the 
middle  station.  This  result  contradicts  the  variability 
imbalance  of  Rao,  which  says  that  the  middle  station  should 
have  a higher  mean  service  time  (case  l.b).  Since  in  Rao's 


-I  c n :r  o c:  Q 3)  :i m<— .-loorm^o 
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Figure  4.1 


Comparison  of  case  lA  and  case  IB 
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formulation  the  variance  of  each  station  is  not  kept  con- 
stant, it  can  be  concluded  that  his  result  is  due  mainly  to 
the  change  in  the  variances  of  service  times  that  is  coupled 
to  the  mean  service  times.  This  conclusion  can  be  explained 
as  follows.  In  Rao's  formulation,  as  the  two  outside 
stations  are  made  faster  their  variances  are  also  made  less. 
Since  the  variance  of  the  middle  station  is  always  zero,  the 
total  variance  is  also  made  less,  and  therefore  one  has  a 
better  production  line.  However,  as  the  two  outside 
stations  become  faster,  the  effect  of  blocking  and  starving 
becomes  more  severe  and  makes  the  throughput  rate  decrease. 
For  this  reason,  the  line  should  be  unbalanced  in  the 
fashion  of  an  inverted  bowl. 

Figure  4.2  shows  contours  of  constant  throughput  rate 
obtained  for  case  2 computed  using  SAS  package.  The  values 
associated  with  the  contours  are  the  relative  throughput 
rate,  which  is  equal  to  one  for  the  balanced  case.  From 
this  figure  we  can  draw  some  conclusions.  First,  when 
y = kept  equal  to  1,  i.e.  ~ ~ ^3  ~ ^ 

are  unbalanced,  the  maximum  improvement  over  the  balanced 
condition  is  0.4%  and  occurs  at  x = As  a 

comparison  Hillier  and  Boling  obtained  an  improvement  of 
about  .52%  at  X = v-^/v2  ~ 1*27  and,  therefore,  y = ^-^1^2  ~ 
1.27.  This  tells  us  that  the  bowl  phenomenon  is  more 
pronounced  when  the  variabilities  are  unbalanced  together 
with  the  mean  values.  We  see  from  the  contour  plot  that  the 
maximum  improvement  is  about  .55%  and  occurs  at  vij^/y2  “ 1.27 
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Figure  4.2.  Result  of  case  2 


Ill 


and  ~ ^•23.  Secondly,  in  order  to  get  the  highest 

throughput  rate  we  should  assign  less  work  (lower  mean 
value)  to  the  more  consistent  (lower  variability)  stations. 
This  trend  is  in  line  with  the  result  of  case  l.a,  but 
contradicts  the  result  of  case  l.b.  Thirdly,  there  are 
points  in  the  contours  at  which  the  throughput  rate  can  be 
increased  by  unbalancing  the  mean  service  times,  but  is 
insensitive  to  unbalancing  of  the  service  time  variabil- 
ities. Then  there  are  other  points  at  which  the  throughput 
rate  can  be  increased  by  unbalancing  the  service  time 
variabilities  but  is  insensitive  to  unbalancing  of  the  mean 
service  times.  In  other  words,  there  are  points  at  which 
the  gradient  of  the  throughput  rate  is  in  the  direction  of 
the  x-axis,  and  there  are  points  at  which  the  gradient  is  in 
the  direction  of  the  y-axis. 

Contours  of  constant  throughput  rate  obtained  for 
case  3 are  given  in  figure  4.3.  This  figure  is  very  similar 
to  figure  4.2  for  case  2 except  that  the  maximum  improvement 
is  .21%  and  occurs  at  x = = 1.09  and  y = ^ ° 2 ~ 

In  figure  4.4  and  4.5  we  show  contours  of  constant 
throughput  rate  for  case  2 and  case  3 applied  to  a paced 
production  line.  We  can  see  from  these  figures  that 
unbalancing  a paced  production  line  in  any  direction 
decreases  the  throughput  rate. 
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Figure  4.3.  Result  of  case  3 
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Figure  4.4.  Result  of  case  2 in  a paced  production  line 
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Figure  4.5.  Result  of  case  3 in  a paced  production  line. 
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SUMMARY 


The  research  presented  in  this  dissertation  was  under- 
taken because  of  the  paucity  of  numerical  methods  that  are 
applicable  to  manufacturing  systems  in  general  and  produc- 
tion lines  in  particular.  The  objective  was  twofold; 
first,  to  develop  an  efficient  solution  procedure  for  the 
throughput  rate  of  a production  line  with  an  arbitrary 
number  of  stations  that  can  handle  a more  general  dis- 
tributions of  service  time,  repair  time  and  time  between 
breakdown  than  the  available  methods;  secondly,  to  use  this 
procedure  to  investigate  the  effect  of  changing  station 
parameter  on  the  line's  throughput  rate. 

Since  in  the  Markov  model  the  distribution  of  service 
time  is  restricted  to  exponential  distribution,  the  funda- 
mental model  used  in  this  dissertation  is  a holding  time 
model.  This  approach  can  be  applied  to  production  lines 
with  an  arbitrary  number  of  stations  and  an  arbitrary 
distribution  of  service  time. 

Following  is  a list  of  the  accomplishments  of  this 
dissertation  in  the  order  of  their  importance. 

1.  We  developed  a solution  procedure  for  the  throughput 
rate  that  is  more  efficient  and  can  handle  a class  of 
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service  time  distributions  that  is  more  general  than 
previous  methods.  We  call  this  class  of  distributions 
special  phase  type  distribution. 

2.  We  obtained  a closed  form  expression  for  the  throughput 
rate  of  a three-station  production  line  that  is  more 
general  than  expressions  heretofore  available. 

3.  We  extended  the  holding  time  model  to  incorporate 
breakdown  and  repair.  This  was  accomplished  through  an 
equivalent  service  time,  called  virtual  service  time, 
that  can  be  used  in  a production  line  without  station 
breakdown.  We  obtained  a general  formula  that  relates 
the  distribution  of  virtual  service  time  to  the  distri- 
butions of  service  time,  repair  time,  and  time  between 
breakdowns . 

Using  these  results,  we  have  been  able  to  do  the  following; 

1.  We  calculated  the  throughput  rate  of  K-station  produc- 
tion lines  without  station  breakdown  for  K = 3,  4, 

. . .,  10  with  general  Erlang  service  time.  The  most 
complex  case  for  which  a solution  was  previously  given 
in  the  literature  is  K = 6 with  exponential  service 
time. 

2.  We  calculated  the  throughput  rate  of  production  lines 
without  station  breakdown  whose  service  time  distribu- 
tions are  in  the  class  of  special  phase  type  that  is 
more  general  than  Erlang  distribution.  Most  analysis  of 
production  lines  previously  carried  by  assuming  special 
Erlang  service  time. 
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3.  We  developed  a method  to  calculate  the  throughput  rates 
of  production  lines  with  station  breakdown  for  the  case 
where  the  distributions  of  service  time,  repair  time  and 
time  between  breakdowns  are  not  exponential.  Previ- 
ously, only  the  case  where  all  three  distributions  are 
exponential  could  be  handled. 

4.  We  produced  an  empirical  formula  for  the  throughput  rate 
of  balanced  production  lines. 

5.  We  developed  a method  of  parametric  analysis  for  the 
throughput  rate  of  production  lines  that  uses  the  mean 
values  and  the  variances  of  service  times  as  the  input. 
Using  this  analysis  we  produced  sets  of  contours  of 
constant  throughput  rate  of  three-station  production 
lines.  These  contours  can  be  used  to  separate  the 
effect  of  mean  value  unbalance  from  that  of  variability 
unbalance  on  the  throughput  rate. 

The  following  open  problems  are  left  as  areas  of  future 
research: 

1 

I j ^ (n)  are  statistically  independent  in  a more  defini- 
tive form. 

2.  Extending  the  current  solution  procedure  in  the  follow- 
ing ways: 

Generalize  the  procedure  so  that  it  can  handle 
service  time  distribution  with  repeated  and/or 
complex  conjugate  poles.  This  case  was  excluded  to 


Solving  the  question  of  whether  or  not  (n-1)  and 
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simplify  the  derivation  of  the  equations  used  in  the 
procedure. 

- Develop  a procedure  that  includes  buffers  in  the 
analysis.  The  current  procedure  treats  buffers  as 
stations  with  zero  service  time.  While  it  is  possi- 
ble to  handle  buffers  this  way,  the  procedure  becomes 
very  inefficient  if  buffers  have  a large  capacity. 
Exploit  the  structure  of  the  resulting  equations  so 
that  they  can  be  solved  more  efficiently.  Prelimi- 
nary investigation  shows  that  the  system  of  equations 
that  need  to  be  solved  has  upper  triangular  struc- 
ture. This  structure  suggests  a better  procedure  to 
solve  this  system  of  equations  than  the  procedure 
currently  used. 

3.  Producing  an  empirical  formula  that  relates  the  through- 
put rate  to  the  mean  and  variance  of  each  station. 
Special  form  of  such  a formula  is  given  in  Chapter  2 for 
the  case  where  all  stations  have  the  same  mean  value  and 


variance. 


APPENDIX  A 


PROBABILITY  DISTRIBUTION  OF  SPECIAL  PHASE  TYPE 

Consider  a Markov  process  with  one  absorbing  state. 

Let  state  i,  i = i,2.  . . . ,m  be  the  transient  states  of 
this  process  and  state  m + 1 be  the  absorbing  state.  Such  a 
process  has  an  (m+1 ) x (m+1 ) transition  rate  matrix  of  the 
following  form; 


(A-1) 


where  the  m x m submatrix  C satisfies  the  following 


properties : 


c . . 

11 


< 0 


S 0 for  i / j 


m 


Also 


Cu  + c = 0 
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where  u is  a vector  with  all  elements  equal  to  1. 

If  this  process  starts  out  in  state  i with  probability 
q^,  i=l,  2,  . . .,m+l,  then  the  distribution  of  time  to 
absorbtion  into  state  m + 1 is  given  as 

[ 1 - q e*“^u  for  t S 0 

F(t)  = < 

( 0 otherwise  (A-2) 

where  q is  a vector  with  q^ , i = 1,  2,  . . . , m as  its  i-th 
element. 

The  probability  distribution  of  phase  type  is  defined 
by  Neuts  [1981]  as  the  probability  of  the  time  to  absorbtion 
in  a Markov  process  with  one  absorbing  state.  Therefore,  in 
general,  this  probability  distribution  can  always  be  written 
in  the  form  (A-2).  The  Laplace-Stielt jes  transform  of  (A-2) 
is 


F*  (S)  = + q (Sl-C)"^  (-C)U 

= + q (sI-C)“^  c (A-3) 

Considering  equation  (A-3) , the  probability  distribution  of 
phase  type  has  the  property  that  its  Laplace-Stielt jes 
transform  is  a rational  function  whose  poles  are  in  the  left 
hand  plane.  Conversely,  a probability  distribution  whose 
Laplace-Stielt jes  transform  is  a rational  function  is  of 
phase  type  that  can  always  be  written  in  the  form  (A-2) . 

This  class  of  probability  distributions  is  a proper  subset 
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of  the  class  of  Coxian  distributions  introduced  by  Cox 
[1955] . 

A subclass  of  the  class  of  phase  type  distributions 
comprises  those  distribution  functions  whose  Laplace- 
Stieltjes  transforms  have  only  real  and  distinct  poles.  We 
call  these  distributions  special  phase  type  distributions. 
If  F(t)  is  a special  phase  type  distribution,  we  can  write 
F(t)  as 


m _h 

1 - E a.  e~°i^  for  t S 0 
i=l  ^ 


F(t)  = 


Otherwise 


(A-4) 


m 

where  b.  >0,0S  E a.  ^1  and  m ^ 1 , and  where  the 
^ i=l  ^ 

parameters  a^  are  restricted  to  values  that  make  F(t)  a 
valid  distribution  function.  The  relation  between  the 
parameters  in  (A-4)  and  (A-2)  is  obtained  by  comparing  the 
Laplace-Stielt jes  transforms  of  both  expressions.  As  a 
result,  we  have 


^m+1 


m 

1 - E 
i=l 


a . 

1 


-b^  is  the  i-th  eigenvalue  of  C 

a.  = lim  (s+b.)  q (sI-C) 

^ s^-b.  ^ 

1 


c 


The  class  of  special  phase  type  distribution  is  quite 
general.  It  includes  general  Erlang  and  hyperexponential 
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distributions  as  special  cases.  The  restriction  to  real  and 
distinct  poles  is  important  for  the  development  of  simple 
numerical  solution  procedures  for  the  problem  discussed  in 
this  study.  Several  properties  of  special  phase  type 
distributions  are  stated  as  theorems  and  lemmas  in  what 
follows . 

Definition:  F„(t)  is  a special  phase  type  distribution  if 


From  this  definition  it  follows  that  X must  be  a nonnegative 
random  variable. 

Theorem  A. 1 i Let  Z = max  [X,Y] . If  F^(t)  and  FY(t)  are 


and  only  if  it  can  be  expressed  in  the  form 


(A-4) . 


special  phase  type  distributions  and  X is 


statistically  independent  of  Y,  then  ^^(t) 
is  also  a special  phase  type  distribution. 


Proof : 


Let  F^it) 


m _h 

1 - Z a.  e”^i^  for  t S 0 
i=l  ^ 


F^it) 


XI  +* 

1 - Z c.  e ^i^  for  t ^ 0 
i=l  ^ 


Then  F^Ct)  = Fjj(t)  Fy(t) 


1 - Z 
i=l 


m 


n 

Z 

i=l 


m 


+ ( I a. 
i=l  ^ 


i=l 
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which  has  a rational  Laplace-Stielt jes 
transform  whose  poles  are  all  distinct,  real 
and  negative. 

Theorem  A. 2;  .Let  Z = max  [X-Y,0].  If  F^(t)  is  a special 
phase  type  distribution  and  Y is  a nonnega- 
tive random  variable  independent  of  X,  then 

F (t)  is  also  a special  phase  type  distri- 

z 

bution. 


Proof;  “ 

F^lt)  = 1 - ; Fy(x)  dF^(x+t)  for  t > 0 


= 1 


Z a.  ( f F (x)b.e“^i^dx)e”^i^ 
i=l  ^ x=0 


= 1 


m 

E 

i=l 


-b . t 

a . g . e i 
1 


which  is  a special  phase  type  distribution. 
Consider  the  following  recursive  relationships  which 
are  given  as  equation  (2-3-1) , (2-3-6)  and  (2-3-7)  in 
Chapter  2 . 


(n)  = 

max 

[H^ 

(n) 

- (n+1) , 0] 

(A- 

■5) 

H . 
3 

(n)  = 

max 

[S. 

(n) 

, Rj^^(n-l)] 

(A- 

■6) 

H . 

(n)  = 

S . (n 

) + 

B . 

(n) 

(A- 

■7) 

3 

3 

3 

= 0 
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where  Sj (n) , jsK  are  random  variables  having  special  phase 
type  distributions.  We  have  the  following  lemmas. 


Lemma  A.l  F , (t)  and  F„  (t)  for  all  je<  are  special  phase 

' ■■  n . 


type  distributions. 


Proof;  Since  (n)  = 0 , F (t)  = u(t)  which  is  a 

^K+1 

special  phase  type  distribution.  Therefore,  from 

equation  (A-6)  for  j=K  and  theorem  A.l,  F (t)  is 

K 

a special  phase  type  distribution.  Using  this 
result  together  with  equation  (A-5)  for  j = K and 
theorem  A. 2,  F (t)  is  also  a special  phase  type 

distribution.  Continuing  in  this  fashion,  we  can 

show  that,  for  all  je<,  F ^(t)  and  Fj^  (t)  are 

Rj  j 

special  phase  type  distributions. 

Lemma  A. 2;  F„  (t)  is  a special  phase  type  distribution. 

D 


From  equations  (A-6)  and  (A-7)  we  have 
Bj  (n)  = max  [Rj_|_^(n-1)  - (n)  , 0] 


Proof; 
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Therefore,  by  theorem  A. 2,  F (t)  is  a special 

D • 


j 


phase  type  distribution. 


Consider  the  recursive  relationships 


(n)  = max  [Sj_^(n)  - (n-1) , 0] 

Rj  (n)  = Hj  (n)  - (n+1) 


(A-8) 


(A. 9) 


which  are  given  as  equation  (2-3-6)  and  (2-3-8)  in 


Chapter  2.  Since  R^ (n)  is  not  a nonnegative  random 


variable,  theorem  A. 2 is  not  applicable  to  equation  (A-8). 
Therefore,  we  cannot  say,  in  general,  that  (t)  will  have 

j 

a special  phase  type  distribution.  It  turns  out  that,  in 
some  cases,  the  Laplace-Stielt jes  transform  of  F^  (t)  has 

repeated  poles.  However,  in  the  approximation  method  we  use 
the  following  relationship 


ij  (n)  = max  [Sj_^(n)  - Rj(n-l),  0] 


(A-10) 


which  is  given  as  equation  (2-6-27)  in  Chapter  2 to 


approximate  the  true  value  of  (n) . We  have  the  following 
lemma. 


Fj  (t)  is  a special  phase  type  distribution. 


Lemma  A . 3 : 
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Proof:  Since  (t)  is  a special  phase  type  distribu- 

tion,  the  proof  follows  from  theorem  A. 2 and 
equation  (A-10) . 

These  three  lemmas  show  that  if  F„  (t)  for  all  je<  are 
special  phase  type  distributions,  then  F ^(t) , (t) , 

F„  (t)  and  F"  (t)  are  special  phase  type  distributions, 

B.  I. 

but  F (t)  is  not.  Since  F~  (t)  is  used  in  the  approxima- 

j j 

tion  method  as  opposed  to  F^  (t)  in  the  exact  method,  these 

lemmas  explain  why  the  approximation  method  is  simpler  and 
more  desirable  than  the  exact  method. 


> 


APPENDIX  B 


A FORMULA  FOR  THE  THROUGHPUT  RATE  OF  A THREE-STATION 

PRODUCTION  LINE 


The  system  of  integral  equations  for  the  holding  time 
model  of  a three-station  production  line  has  a special 
structure.  For  the  case  when  the  service  time  distributions 
of  the  two  outside  stations  are  of  special  phase  type  this 
system  of  integral  equations  can  be  transformed  into  a 
system  of  linear  algebraic  equations  with  the  number  of 
equations  equal  to  the  number  of  stages  of  the  service  time 
distribution  of  the  last  station.  Furthermore,  closed  form 
expressions  for  the  throughput  rate  can  also  be  obtained  for 
special  service  time  distributions  of  the  last  station. 

When  K = 3,  equations  (2-4-4)  through  (2-4-6)  reduce 
to 

00 

/ F„  (X)  F _^(x)  dF„  (x+t)  for  t > 0 
x=0  ^2  R3  1 

otherwise  (B-1) 

00 

/ F (X)  dF  (x+t)  for  t > 0 

x=0  ^2  ^3 

otherwise  (B-2) 
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Let  the  service  time  distributions  of  station  1 and  3 
be  of  special  phase  type,  i.e. 


m -a  . t 

F„  (t)  = 1 - Z B.e  ^ 
^1  i=l  ^ 


for  t S 0 (B-3) 


and 


(t) 

^3 


n 

1 - Z 

j = l 


for  t > 0 (B-4) 


where  m and  n are  the  number  of  stages  of  these  distribu- 
tions . 


Substituting  (B-4)  into  (B-2)  yields 


F , (t)  = 1 - Z if  F (X)  6.y.e  dx)  e 
R3  j = l 0-^2  J ^ 


n -Y -t 

= 1 - Z a.e  for  t S 0 (B-5) 

j = l ^ 


where 


00  ~Y  • X 

a.  — f F^  (x)  Y’^^®  ^ forj=l/2,  . • .,n 

3 0 2 ^ ^ 


(B-6) 
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Similarly,  substituting  (B-3)  into  (B-1)  yields 


m » -a . X 

(t)  = 1 - E ( / F„  (X)  F . (X)  S.a.e  ^ dx)e 
^2  i=l  0 ^2  R3  ^ ^ 


m -a  . t 

= 1 - E b.e  ^ for  t £ 0 

i=l  ^ 


where 


® -a . X 

b.  = / F„  (x)  F (x)  g.a.e  ^ dx 

^ 0 ^2  R3  ^ ^ 


for  i = 1 , 2 , . . . , m 


Substituting  (B-7)  into  (B-6)  yields 


a . 


3 


6 . 
3 


m 

E 

i=l 


b . Y • 5 . 
I'l  3 


/ (otj_  + Yj) 


for  j = 1 , 2 , . . . , m 


and  substituting  (B-5)  into  (B-8)  yields 


00 

b.  = / F„  (X)  (1 

^ 0 ^2 


n ~Y • X -a . X 

E a.e  ^ ) 6-a.e  ^dx 

j = l 3 


= 6.cx.F3^(a.)  - 


_E^  aj6.a.Fg^(a.  + Yj) 


(B-7) 


(B-8) 


(B-9) 


for  i 


1,  2, 


• f 


m 


(B-IO) 
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where  (s)  is  the  Laplace  transform  of  F (t) . 

S2  ^2 

To  put  equations  (B-9)  and  (B-10)  in  a more  compact 

form,  we  define  the  column  vectors 


a = 

[a^,  a^^]^ 

b = 

^2 ' * ' 

c = 

[Ci»  ^2'  • ' ■'  ‘^m^ 

where  Cj^ 

= B.a.  F„  (a.).  We  further  define  the  matrices 
1 1 S2  1 

P = 

[p.  • ] 

j mxn 

Q = 

[q  • ■ ] 

‘■^1]  nxm 

where 

Pij 

= 6.0.  Fg^(c.  + T.) 

= / l»i  + Tj) 

Equations  (B-9)  and  (B-10)  becomes 


a = 

6 - Qb  (B-11) 

b = 

c — Pa  (B— 12 ) 

where  6 is  a column  vector  whose  j-th  element  is  6^. 
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Eliminating  b from  these  two  equations  yields  a system 
of  equations  defining  the  vector  a 


a = 6 - Q (c  - Pa) 


or 


a = (I  - QP)“^  (6  - Qc) 


(B-13) 


Given  the  parameters  a^,  y j / <5^  and  the  Laplace  trans- 

form of  the  service  time  distribution  of  the  middle  station, 
the  value  of  a is  calculated  with  equation  (B-13).  Once 
this  value  has  been  found,  the  mean  interdemand  time  E[H^] 
is  calculated  by  substituting  (B-5)  into  (2-5-5)  and  then 
(2-5-3);  namely 


E[H,]  = / (1  - (t)  F (t)  F (t))dt 

^0  ^1  ^2  R3 


00  n -Y't 

= / (1  - F„  (t)  F„  (t)  (1  - Z a.e  ))dt 

0 ^1  ^2  j=l  ^ 


Substituting  (B-3)  for  F (t)  in  the  above  equation  yields 

^1 
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m 


-a  . t 
1 


E[HJ  = f {1  - F„  (t))dt  + S g.  / F (t)e  dt 
^0  ^2  i=l  ^ 0 ^2 


n ® m -a.t  “Y-;t 

+ Z a.  ( / F„  (t)  {1  - Z g.e  ^ )e  ^ dt) 

j=l  ^ 0 ^2  i=l  ^ 


m 

= u-j  + 2 B-F„  (a.)  + 

i=l  1 ^ 


n 

2 a.(F  (Y.) 
j=l  ^ ^2  ^ 


m 

Z 

i=l 


''j” 


(B-14) 


where  U2  is  the  mean  service  times  of  the  middle  station. 

From  the  above  derivation  we  can  see  that  formula 
(B-14)  can  be  used  for  arbitrary  service  time  distribution 
of  the  middle  station  as  long  as  the  Laplace  transform  of 
this  distribution  function  is  known.  For  example,  if  the 
middle  station  has  a deterministic  service  time  equal  to  ^2/ 

then  F (s)  = e /s. 

^2 


Equation  (B-14)  is  not  a closed  form  expression  since 
in  order  to  find  the  value  of  a we  must  invert  the  nxn 
matrix  (I-QP) . In  general  this  inversion  must  be  done 
numerically.  However,  for  some  special  cases,  the  value  of 
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a is  readily  available  from  equation  (B-13) . These  special 
cases  are  given  below 

Case  1 ; m = n = 1 

This  case  means  that  the  two  outside  stations  have 
exponential  service  times.  Therefore 

= 6j  = 1 ; = a and  Yj  = Y 

Equation  (B-14)  reduces  to 


where  a is  obtained  from  equation  (B-13)  as  follows 

a = (1  - qc)  / (1  - qp) 
c = a F„  (a) 

p = a F-  (a  + y) 
q = Y / (a  + y) 

This  closed  form  expression  is  more  general  than  the 


E[H^]  = ^2  + Fg  (a)  + a (F 


F„  (a  + y) ) 
^2 


(B-15) 


ones  given  by  Hunt  [1956],  Rao  [1976]  or  Muth  [1984]  in  the 
sense  that  the  middle  station  can  have  arbitrary  service 
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time  distribution.  While  Rao  [1976]  derived  two  different 
formulas  for  the  middle  station  having  deterministic  and 
uniform  service  times,  and  Muth  [1984]  derived  a formula  for 
the  middle  station  having  exponential  service  time,  these 
three  cases  are  represented  here  by  the  general  formula 
given  by  (B-15) . 

Case  2;  n = 1 

In  this  case  only  the  last  station  has  exponential 
service  time.  Therefore 

= 1 and  Yj  = Y 

Equation  (B-14)  reduces  to 

m 

E[H.]  = y.  + I e.  F (a.)  + 

^ ^ i=l  ^ ^ 

m 

a (F„  (y)  - 2 e.  F„  (a.  + y))  (B-16) 

^2  i=l  ^ ^2  ^ 


where 
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Case 


case 


m 


m 


a = (1  - E 2 

i=l  i=l 


' “i®i 


Pi  = “I®!  Psj'”!  + I"* 


= Y / («£  + y) 


3:  n = 2 


Without  going  into  details  of  the  derivations,  for  this 
we  obtain  from  equation  (B-13) 


m 


m 


*1  = I <1  - ‘521P12'  - .fj  ‘Jll=i> 


m 


m 


+ I 9iiPi2)  <«2  - ■321=1* 


m 


m 


*2  = 


= A 13  - 3llPll*  <*2  - 32iCi> 


m 


m 


^ A '^2iPil^  "^li^i^ 

1=1  1=1 


where 


A = (1  - 


m 

Z 

i=l 


) (1  - 


m 

Z 

i=l 


'J2iPl2' 


m 

- ( 2 
i=l 


'3llPi2> 


m 
( E 
i=l 


^2kPi2 


) 


Pij  *^ij  given  previously.  The  mean 

interdemand  time  is  calculated  using  equation  (B-14)  with 


and  as  given  above. 


APPENDIX  C 


THE  DERIVATION  OF  F_  (t)  AND  F (t) 

I R. 


We  derive  the  formula  for  the  distribution  of  a random 
variable  that  is  a function  of  the  difference  of  two  inde- 
pendent random  variables.  Let  X and  Y be  two  independent 
random  variables  with  distribution  functions  F^^(t)  and 
F^(t) . The  distribution  function  of  X - Y can  be  derived 
using  the  law  of  total  probability.  By  conditioning  on 
X = X,  we  obtain 


F^_^{t)  = P[X-Y  ^ t] 


X = -00 


/ P[x  - Y S t]  d Fj^(x) 


; (1  - FY(x-t)  ) d Fj^(x) 

X=-oo 


1 - f Fylx-t)  d Fj^(x) 

X = -00 


(C-1) 
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If  X is  a nonnegative  random  variable,  then 
(C-1)  reduces  to 


F„  ^(t)  = 1 - / F (x-t)  d F (X) 

x=0  ^ 


By  letting  x-t  = y,  equation  (C-2)  becomes 


Fjj_y(t)  = 1 - ; Fy(Y)  d Fj^(y  + t) 


If  Y is  also  a nonnegative  random  variable,  then 
(C-2)  reduces  to 


1 - / F (x-t)  d F^(x)  for 

x=t 


Fx_y(t)  = 


1 - f Fy(x-t)  d Fjj(x)  for 

x=0 


and  equation  (C-3)  reduces  to 


1 - 


f F^(y)  d F^(y  + t) 
y=0 


1 - 


/ F (y)  d F^iy  + t) 
y=-t 


for 


equation 


(C-2) 


(C-3) 


equation 


t ^ 0 


t < 0 (C-4) 


for  t > 0 


t < 0 


(C-5) 
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Let  the  random  variable  Z be  defined  as  Z = max [X-Y , 0] 
where  X is  a nonnegative  random  variable,  then  the  distribu- 
tion function  of  Z is  obtained  as 


F^(t) 


/ Fyly)  d Fjj(y  + t)  for  t ^ 0 


for  t < 0 


(C-6) 


Consider  the  random  variable  Rj  that  is  computed  using 
relationship  (n)  = (n)  - Ij_j^(n+1)  given  as  equation 

(2-3-6)  in  Chapter  2,  where  (n)  and  Ij_^(n+1)  are  both 
nonnegative  random  variables.  From  equations  (C-4)  and 
(C-5)  we  obtain 


F (t) 


^ 1 - f F^  (x)  d F„  (X  + t)  for  t > 0 
x=0  ^j-1  “j 


^1  - F (x-t)  d Fjj  (X)  for  t < 0 

x=0  j-1  j 

(C-7) 


Since  Rj(n)  = max[Rj(n),  0]  and  R^  (n)  = min[Rj(n),  0],  we 


have 
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1 - / (x)  d F„  (x  + t)  for  t S 0 

x=0  ^j-1  “j 


F (t)  = 

rT 

D 


for  t < 0 

(C-8) 


and 


( 1 - / F^  (x-t)  d F„  (x)  for  t ^ 0 

x=0  D-1  D 


_(t)  = / 


for  t > 0 

(C-9) 


Similarly,  the  distribution  function  of  Ij  is  obtained  from 
the  relationship  Ij (n)  = max[Sj_^(n)  - (n-1) , 0]  given  as 

equation  (2-3-8)  in  Chapter  2,  where  Sj_^  is  a nonnegative 
random  variable.  From  equation  (C-6)  we  have 


Fj  (t)  = 

j 


1 - / F_  (X)  d F„  (X  + t)  for  t S 0 

x=-t  ^j-1 


for  t < 0 
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or 


/I  - / F (x)dF  (x+t)  - 


x=-t  Rj 


/ F (x)dF  (x+t) 
x=0  Rj  S._^  . 


for  t ^ 0 


Fj  (t) 


V 


for  t < 0 


(C-10) 
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